GESTIÓ DELS ESPAIS VERDS PER LA CONSERVACIÓ DE POL·LINITZADORS:
Carregar llibreries
#install.packages("vegan")
library(openxlsx)
library(readxl)
library(dplyr)
library(vegan)
library(MASS)
library(reshape2)
library(reshape)
library(ggplot2)
library(eulerr)
library(purrr)
library(stringr)
library(ggbiplot)
library(indicspecies)
library(purrr)
library(lme4)
library(Matrix)
library(lmerTest)Carregar la base de dades original (corregida) i eliminar les dades de Mordèl·lids i de les zones “Llig”, “Mura”, i “TV-7041”
# Obtenir el directori de l'arxiu actual
directorio <- getwd()
dades_finals <- read_excel(file.path(directorio,"Carreteres_ITS+altres.xlsx"))
dades_finals <- dades_finals %>%
filter(Familia != "Mordellidae" & !Zona %in% c("Llig", "Mura", "TV-7041"))
head(dades_finals)## # A tibble: 6 × 25
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 28 Campus Juliol 7 2021 4 NS <NA>
## # ℹ 17 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>ANÀLISI DESCRIPTIVA
# Definir nova variable: "Mostreig", segons Zona, Mes i Any
dades_finals <- dades_finals %>%
mutate(Mostreig = paste(Zona, Mes, Any, sep = "-"))Riquesa
Riquesa total
## # A tibble: 1 × 1
## riquesa_total
## <int>
## 1 137Riquesa per cada mostreig i tractament
## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.## # A tibble: 11 × 3
## # Groups: Mostreig [6]
## Mostreig Tractament riquesa
## <chr> <chr> <int>
## 1 Campus-Juliol-2021 NS 20
## 2 Campus-Juliol-2021 S 21
## 3 Campus-Jun-2021 NS 45
## 4 Campus-Jun-2021 S 34
## 5 Campus-juliol-2020 NS 39
## 6 Campus-juliol-2020 S 24
## 7 Campus-maig-2020 NS 45
## 8 Franqueses-Maig-2021 NS 20
## 9 Franqueses-Maig-2021 S 11
## 10 Sabadell-Juny-2021 NS 16
## 11 Sabadell-Juny-2021 S 12Riquesa per cada zona i tractament
riquesa_zones <- dades_finals %>%
group_by(Zona, Tractament) %>%
summarise(riquesa = n_distinct(ID))## `summarise()` has grouped output by 'Zona'. You can override using the
## `.groups` argument.## # A tibble: 6 × 3
## # Groups: Zona [3]
## Zona Tractament riquesa
## <chr> <chr> <int>
## 1 Campus NS 98
## 2 Campus S 59
## 3 Franqueses NS 20
## 4 Franqueses S 11
## 5 Sabadell NS 16
## 6 Sabadell S 12Riquesa per grups
# Abelles
riquesa_abelles<- dades_finals %>%
filter(Ordre == 'HymenopteraAbella') %>%
summarise(riquesa_abelles=n_distinct(ID))
riquesa_abelles## # A tibble: 1 × 1
## riquesa_abelles
## <int>
## 1 80# Coleòpters
riquesa_coleoptera<- dades_finals %>%
filter(Ordre == 'Coleoptera') %>%
summarise(riquesa_coleoptera=n_distinct(ID))
riquesa_coleoptera## # A tibble: 1 × 1
## riquesa_coleoptera
## <int>
## 1 33# Sírfids
riquesa_sirfids<- dades_finals %>%
filter(Familia == 'Syrphidae') %>%
summarise(riquesa_sirfids=n_distinct(ID))
riquesa_sirfids## # A tibble: 1 × 1
## riquesa_sirfids
## <int>
## 1 9Abundància
Abundància total
## # A tibble: 1 × 1
## Numero_total
## <int>
## 1 1111Abundància per cada zona i tractament
## `summarise()` has grouped output by 'Zona'. You can override using the
## `.groups` argument.## # A tibble: 6 × 3
## # Groups: Zona [3]
## Zona Tractament Numero_total
## <chr> <chr> <int>
## 1 Campus NS 649
## 2 Campus S 323
## 3 Franqueses NS 45
## 4 Franqueses S 14
## 5 Sabadell NS 29
## 6 Sabadell S 51Abundància de cada gènere total
## # A tibble: 6 × 2
## Genere Numero_total
## <chr> <int>
## 1 Acmaeodera 6
## 2 Amegilla 1
## 3 Andrena 6
## 4 Anthaxia 29
## 5 Anthidium 10
## 6 Apis 11Abundància de cada ID per zona i tractament
abundaciazona_ID <- dades_finals %>%
group_by(Zona, Tractament, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Zona', 'Tractament'. You can override
## using the `.groups` argument.## # A tibble: 6 × 4
## # Groups: Zona, Tractament [1]
## Zona Tractament ID Numero_total
## <chr> <chr> <chr> <int>
## 1 Campus NS Acmaeodera cylindrica 5
## 2 Campus NS Amegilla sp. 1
## 3 Campus NS Anthaxia (Anthaxia) 1
## 4 Campus NS Anthaxia (Melanthaxia) 5
## 5 Campus NS Anthaxia hungarica 10
## 6 Campus NS Anthaxia sp. 7Abundancia de cada ID per cada mostreig i tractament
abundacia_ID <- dades_finals %>%
group_by(Mostreig, Tractament, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.## # A tibble: 6 × 4
## # Groups: Mostreig, Tractament [1]
## Mostreig Tractament ID Numero_total
## <chr> <chr> <chr> <int>
## 1 Campus-Juliol-2021 NS Amegilla sp. 1
## 2 Campus-Juliol-2021 NS Anthaxia sp. 1
## 3 Campus-Juliol-2021 NS Bombus pascuorum 1
## 4 Campus-Juliol-2021 NS Ceratina cucurbitina 1
## 5 Campus-Juliol-2021 NS Ceratina parvula 7
## 6 Campus-Juliol-2021 NS Chrysididae 1PLOTS ABUNDÀNCIES PER GÈNERE
Campus
ab_abellescampus <- dades_finals %>%
filter(Zona == 'Campus') %>%
filter(Ordre == 'HymenopteraAbella') %>%
group_by(Genere, Tractament) %>%
summarise(Numero_total = n(), .groups = "drop")
p <-ggplot(ab_abellescampus, aes(y = reorder(Genere, Numero_total), x = Numero_total, fill = Tractament)) +
geom_col(position = position_dodge(preserve = "single"),width = 0.75) +
geom_label(aes(label = Numero_total,x= -5,y=ifelse(Tractament=="NS",as.numeric(reorder(Genere, Numero_total))-0.25,as.numeric(reorder(Genere, Numero_total))+0.2)),
position = position_dodge(width = 0.3),
show.legend = FALSE,
size = 3.5,
fill = ifelse(ab_abellescampus$Tractament=="NS","#81b29a","#e07a5f")) +
scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) + # Ajusta los nombres de los tratamientos
guides(fill = guide_legend(nrow = 2, byrow = TRUE)) +
theme_bw() +
theme(legend.position = "right",
axis.text.y = element_text(size = 15, face = "italic"),
axis.text.x = element_text(size = 15),
axis.title.x = element_text(size = 15),
axis.title.y = element_text(size = 15, face = "italic")
) +
expand_limits(x = -6) +
labs(x = "Abundància",
y = "")
p
Sabadell
ab_abellessabadell <- dades_finals %>%
filter(Zona == 'Sabadell') %>%
filter(Ordre == 'HymenopteraAbella') %>%
group_by(Genere, Tractament) %>%
summarise(Numero_total = n(), .groups = "drop")
q <-ggplot(ab_abellessabadell, aes(y = reorder(Genere, Numero_total), x = Numero_total, fill = Tractament)) +
geom_col(position = position_dodge(preserve = "single"),width = 0.75) +
geom_label(aes(label = Numero_total,x= 0,y=ifelse(Tractament=="NS",as.numeric(reorder(Genere, Numero_total))-0.25,as.numeric(reorder(Genere, Numero_total))+0.2)),
position = position_dodge(width = 0.3),
show.legend = FALSE,
size = 3.5,
fill = ifelse(ab_abellessabadell$Tractament=="NS","#81b29a","#e07a5f")) +
scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) + # Ajusta los nombres de los tratamientos
guides(fill = guide_legend(nrow = 2, byrow = TRUE)) +
theme_bw() +
theme(legend.position = "right",
axis.text.y = element_text(size = 15, face = "italic"),
axis.text.x = element_text(size = 15),
axis.title.x = element_text(size = 15),
axis.title.y = element_text(size = 15, face = "italic")
) +
labs(x = "Abundància",
y = "")
q
Franqueses
ab_abellesfranqueses <- dades_finals %>%
filter(Zona == 'Franqueses') %>%
filter(Ordre == 'HymenopteraAbella') %>%
group_by(Genere, Tractament) %>%
summarise(Numero_total = n(), .groups = "drop")
r <-ggplot(ab_abellesfranqueses, aes(y = reorder(Genere, Numero_total), x = Numero_total, fill = Tractament)) +
geom_col(position = position_dodge(preserve = "single"),width = 0.75) +
geom_label(aes(label = Numero_total,x= 0,y=ifelse(Tractament=="NS",as.numeric(reorder(Genere, Numero_total))-0.25,as.numeric(reorder(Genere, Numero_total))+0.2)),
position = position_dodge(width = 0.3),
show.legend = FALSE,
size = 3.5,
fill = ifelse(ab_abellesfranqueses$Tractament=="NS","#81b29a","#e07a5f")) +
scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) + # Ajusta los nombres de los tratamientos
guides(fill = guide_legend(nrow = 2, byrow = TRUE)) +
theme_bw() +
theme(legend.position = "right",
axis.text.y = element_text(size = 15, face = "italic"),
axis.text.x = element_text(size = 15),
axis.title.x = element_text(size = 15),
axis.title.y = element_text(size = 15, face = "italic")
) +
labs(x = "Abundància",
y = "")
r
DIVERSITAT DE SHANNON I CORBES D’ACUMULACIÓ D’ESPÈCIES
Diversitat de Shannon per Zona i Tractament
# Preparar dades: Crear una taula de l'abundància de cada ID (columnes) segons la zona i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux2 <- cast(abundaciazona_ID, Zona + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux2 <- as.data.frame(taula_aux2)
taula_aux2[is.na(taula_aux2)] <- 0# Agrupar Zona amb Tractament (en una nova variable Agrupacio2)
data_shannon_ZT <- taula_aux2 %>%
mutate(Agrupacio2 = paste(Zona, Tractament, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Zona, -Tractament)# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio2 i establir-la com a identificador de les files
rownames(data_shannon_ZT) <- data_shannon_ZT$Agrupacio2
data_shannon_ZT <- data_shannon_ZT %>% dplyr::select(-Agrupacio2) %>% mutate_all(as.numeric)# Càlcul de la diversitat de Shannon per cada Mostreig i Tractament
shannondiv2 <- diversity(data_shannon_ZT)
print(shannondiv2)## Campus_NS Campus_S Franqueses_NS Franqueses_S Sabadell_NS
## 3.612943 2.870658 2.627033 2.341994 2.382088
## Sabadell_S
## 2.003690Diversitat de Shannon per Mostreig i Tractament
# Preparar dades: Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux <- cast(abundacia_ID, Mostreig + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux <- as.data.frame(taula_aux)
taula_aux[is.na(taula_aux)] <- 0# Agrupar Mostreig amb Tractament (en una nova variable Agrupacio)
data_shannon_MT <- taula_aux %>%
mutate(Agrupacio = paste(Mostreig, Tractament, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament)# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_MT) <- data_shannon_MT$Agrupacio
data_shannon_MT <- data_shannon_MT %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)# Càlcul de la diversitat de Shannon per cada Mostreig i Tractament
shannondiv <- diversity(data_shannon_MT)
print(shannondiv)## Campus-juliol-2020_NS Campus-juliol-2020_S Campus-Juliol-2021_NS
## 3.011325 2.404801 2.593731
## Campus-Juliol-2021_S Campus-Jun-2021_NS Campus-Jun-2021_S
## 1.933200 3.161179 2.856531
## Campus-maig-2020_NS Franqueses-Maig-2021_NS Franqueses-Maig-2021_S
## 2.957046 2.627033 2.341994
## Sabadell-Juny-2021_NS Sabadell-Juny-2021_S
## 2.382088 2.003690Corba d’acumulació d’espècies per mostreig (11 mostrejos)
## Warning in cor(x > 0): the standard deviation is zero
Extrapolació de dades (quantes espècies més trobaríem si féssim més mostrejos)
## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 137 236.8148 31.40658 207 26.07681 250.1182 167.3704 12.88879 11
En aquest mostreig s’han trobat 137 espècies (riquesa), però s’estima que hi ha altres espècies “unseen” (no detectades). Segons diferents estimacions, la riquesa total tenint en compte aquestes espècies no detectades seria 236 (chao), 207 (jackknife primer ordre), 250 (jackknife segon ordre), i 167 (bootstrap).
Diversitat de Shannon per Mostreig, Tractament i Trampa
# Preparar dades:
#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa
abundacia_trampa <- dades_finals %>%
group_by(Mostreig, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.## # A tibble: 564 × 5
## # Groups: Mostreig, Tractament, Trampa [101]
## Mostreig Tractament Trampa ID Numero_total
## <chr> <chr> <dbl> <chr> <int>
## 1 Campus-Juliol-2021 NS 2 Lasioglossum glabriusculum 1
## 2 Campus-Juliol-2021 NS 4 Curculionidae 1
## 3 Campus-Juliol-2021 NS 4 Lasioglossum malachurum 1
## 4 Campus-Juliol-2021 NS 4 Pompilidae 1
## 5 Campus-Juliol-2021 NS 4 Pyronia cecilia 1
## 6 Campus-Juliol-2021 NS 6 Lasioglossum glabriusculum 1
## 7 Campus-Juliol-2021 NS 6 Lasioglossum politum 3
## 8 Campus-Juliol-2021 NS 6 Pompilidae 1
## 9 Campus-Juliol-2021 NS 8 Chrysididae 1
## 10 Campus-Juliol-2021 NS 8 Hylaeus sp. 1
## # ℹ 554 more rows#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux <- cast(abundacia_trampa, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_aux <- as.data.frame(taula_aux)
taula_aux[is.na(taula_aux)] <- 0# Agrupar Mostreig amb Tractament i Trampa
data_shannon_MTTR <- taula_aux %>%
mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament, -Trampa)# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_MTTR) <- data_shannon_MTTR$Agrupacio
data_shannon_MTTR <- data_shannon_MTTR %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondiv <- diversity(data_shannon_MTTR)
print(shannondiv)## Campus-juliol-2020_NS_2 Campus-juliol-2020_NS_3 Campus-juliol-2020_NS_4
## 1.7917595 1.8310205 2.0449312
## Campus-juliol-2020_NS_6 Campus-juliol-2020_NS_8 Campus-juliol-2020_NS_10
## 1.7328680 1.8310205 1.0735428
## Campus-juliol-2020_NS_11 Campus-juliol-2020_NS_12 Campus-juliol-2020_NS_13
## 0.0000000 1.7478681 1.6434177
## Campus-juliol-2020_NS_14 Campus-juliol-2020_NS_15 Campus-juliol-2020_NS_16
## 1.1240357 1.5498260 1.3862944
## Campus-juliol-2020_NS_20 Campus-juliol-2020_S_1 Campus-juliol-2020_S_5
## 2.2718685 0.6931472 1.9072840
## Campus-juliol-2020_S_7 Campus-juliol-2020_S_9 Campus-juliol-2020_S_17
## 0.6924212 0.0000000 1.0986123
## Campus-juliol-2020_S_18 Campus-juliol-2020_S_19 Campus-juliol-2020_S_21
## 1.6094379 1.5607104 1.1189689
## Campus-juliol-2020_S_22 Campus-Juliol-2021_NS_2 Campus-Juliol-2021_NS_4
## 1.0751393 0.0000000 1.3862944
## Campus-Juliol-2021_NS_6 Campus-Juliol-2021_NS_8 Campus-Juliol-2021_NS_10
## 0.9502705 1.4648164 0.6931472
## Campus-Juliol-2021_NS_12 Campus-Juliol-2021_NS_14 Campus-Juliol-2021_NS_16
## 1.3296613 1.3321790 0.0000000
## Campus-Juliol-2021_NS_18 Campus-Juliol-2021_NS_20 Campus-Juliol-2021_S_1
## 1.3862944 1.0986123 1.0986123
## Campus-Juliol-2021_S_3 Campus-Juliol-2021_S_5 Campus-Juliol-2021_S_7
## 0.6365142 2.0963526 1.5000764
## Campus-Juliol-2021_S_9 Campus-Juliol-2021_S_11 Campus-Juliol-2021_S_13
## 1.1973401 1.4388830 0.0000000
## Campus-Juliol-2021_S_15 Campus-Juliol-2021_S_17 Campus-Jun-2021_NS_2
## 0.6931472 0.5091373 2.3345491
## Campus-Jun-2021_NS_4 Campus-Jun-2021_NS_6 Campus-Jun-2021_NS_8
## 2.0854684 1.1537419 2.6197294
## Campus-Jun-2021_NS_10 Campus-Jun-2021_NS_12 Campus-Jun-2021_NS_14
## 1.8310205 2.5575077 1.7623137
## Campus-Jun-2021_NS_16 Campus-Jun-2021_NS_18 Campus-Jun-2021_NS_20
## 1.3593685 1.9061547 0.9502705
## Campus-Jun-2021_S_1 Campus-Jun-2021_S_3 Campus-Jun-2021_S_5
## 2.2718685 2.0028830 1.7351265
## Campus-Jun-2021_S_7 Campus-Jun-2021_S_9 Campus-Jun-2021_S_11
## 0.9502705 2.1383972 1.5595812
## Campus-Jun-2021_S_13 Campus-Jun-2021_S_15 Campus-Jun-2021_S_17
## 1.8576598 0.7356219 2.1639557
## Campus-Jun-2021_S_19 Campus-maig-2020_NS_1 Campus-maig-2020_NS_2
## 1.9722470 0.8392696 1.3642812
## Campus-maig-2020_NS_3 Campus-maig-2020_NS_4 Campus-maig-2020_NS_5
## 1.6957425 2.5521719 1.5171064
## Campus-maig-2020_NS_6 Campus-maig-2020_NS_7 Campus-maig-2020_NS_8
## 2.2614951 1.2424533 1.7480673
## Campus-maig-2020_NS_9 Campus-maig-2020_NS_10 Campus-maig-2020_NS_11
## 1.6769878 1.9722470 1.7480673
## Campus-maig-2020_NS_12 Campus-maig-2020_NS_13 Campus-maig-2020_NS_14
## 1.0986123 1.5481155 1.8891592
## Campus-maig-2020_NS_15 Campus-maig-2020_NS_16 Campus-maig-2020_NS_17
## 1.8392967 1.0356097 1.3296613
## Campus-maig-2020_NS_18 Campus-maig-2020_NS_19 Campus-maig-2020_NS_20
## 1.5810938 1.3535914 0.9002561
## Franqueses-Maig-2021_NS_1 Franqueses-Maig-2021_NS_2 Franqueses-Maig-2021_NS_3
## 1.0789922 0.8675632 0.6931472
## Franqueses-Maig-2021_NS_4 Franqueses-Maig-2021_NS_5 Franqueses-Maig-2021_S_1
## 2.2739657 0.7549968 1.3321790
## Franqueses-Maig-2021_S_2 Franqueses-Maig-2021_S_4 Sabadell-Juny-2021_NS_1
## 1.0986123 1.5607104 1.7478681
## Sabadell-Juny-2021_NS_4 Sabadell-Juny-2021_NS_7 Sabadell-Juny-2021_NS_8
## 0.9502705 1.7917595 1.4750763
## Sabadell-Juny-2021_NS_11 Sabadell-Juny-2021_S_2 Sabadell-Juny-2021_S_3
## 1.3862944 0.6365142 0.5623351
## Sabadell-Juny-2021_S_5 Sabadell-Juny-2021_S_6 Sabadell-Juny-2021_S_9
## 1.4995094 1.0986123 0.6931472
## Sabadell-Juny-2021_S_10 Sabadell-Juny-2021_S_12
## 0.5004024 1.2852930Corba d’acumulació d’espècies segons mostreig, tractament i trampa
Extrapolació de dades (quantes espècies més trobaríem si féssim més mostrejos)
## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 137 224.1493 29.51259 201.3564 11.27322 241.7798 164.5183 5.987207 101
Amb aquest mostreig s’han trobat 137 espècies (riquesa), però s’estima que hi ha altres espècies “unseen” (no detectades). Segons diferents estimacions, la riquesa total tenint en compte aquestes espècies no detectades seria 224 (chao), 201 (jackknife primer ordre), 241 (jackknife segon ordre), i 164 (bootstrap).
Corba i extrapolació pel Campus:
#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades pel Campus
abundacia_trampa_campus <- dades_finals %>%
filter(Zona == 'Campus') %>%
group_by(Mostreig, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxcampus <- cast(abundacia_trampa_campus, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxcampus <- as.data.frame(taula_auxcampus)
taula_auxcampus[is.na(taula_auxcampus)] <- 0
# Agrupar Mostreig amb Tractament i Trampa
data_shannon_campus <- taula_auxcampus %>%
mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament, -Trampa)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_campus) <- data_shannon_campus$Agrupacio
data_shannon_campus <- data_shannon_campus %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivcampus <- diversity(data_shannon_campus)
print(shannondivcampus)## Campus-juliol-2020_NS_2 Campus-juliol-2020_NS_3 Campus-juliol-2020_NS_4
## 1.7917595 1.8310205 2.0449312
## Campus-juliol-2020_NS_6 Campus-juliol-2020_NS_8 Campus-juliol-2020_NS_10
## 1.7328680 1.8310205 1.0735428
## Campus-juliol-2020_NS_11 Campus-juliol-2020_NS_12 Campus-juliol-2020_NS_13
## 0.0000000 1.7478681 1.6434177
## Campus-juliol-2020_NS_14 Campus-juliol-2020_NS_15 Campus-juliol-2020_NS_16
## 1.1240357 1.5498260 1.3862944
## Campus-juliol-2020_NS_20 Campus-juliol-2020_S_1 Campus-juliol-2020_S_5
## 2.2718685 0.6931472 1.9072840
## Campus-juliol-2020_S_7 Campus-juliol-2020_S_9 Campus-juliol-2020_S_17
## 0.6924212 0.0000000 1.0986123
## Campus-juliol-2020_S_18 Campus-juliol-2020_S_19 Campus-juliol-2020_S_21
## 1.6094379 1.5607104 1.1189689
## Campus-juliol-2020_S_22 Campus-Juliol-2021_NS_2 Campus-Juliol-2021_NS_4
## 1.0751393 0.0000000 1.3862944
## Campus-Juliol-2021_NS_6 Campus-Juliol-2021_NS_8 Campus-Juliol-2021_NS_10
## 0.9502705 1.4648164 0.6931472
## Campus-Juliol-2021_NS_12 Campus-Juliol-2021_NS_14 Campus-Juliol-2021_NS_16
## 1.3296613 1.3321790 0.0000000
## Campus-Juliol-2021_NS_18 Campus-Juliol-2021_NS_20 Campus-Juliol-2021_S_1
## 1.3862944 1.0986123 1.0986123
## Campus-Juliol-2021_S_3 Campus-Juliol-2021_S_5 Campus-Juliol-2021_S_7
## 0.6365142 2.0963526 1.5000764
## Campus-Juliol-2021_S_9 Campus-Juliol-2021_S_11 Campus-Juliol-2021_S_13
## 1.1973401 1.4388830 0.0000000
## Campus-Juliol-2021_S_15 Campus-Juliol-2021_S_17 Campus-Jun-2021_NS_2
## 0.6931472 0.5091373 2.3345491
## Campus-Jun-2021_NS_4 Campus-Jun-2021_NS_6 Campus-Jun-2021_NS_8
## 2.0854684 1.1537419 2.6197294
## Campus-Jun-2021_NS_10 Campus-Jun-2021_NS_12 Campus-Jun-2021_NS_14
## 1.8310205 2.5575077 1.7623137
## Campus-Jun-2021_NS_16 Campus-Jun-2021_NS_18 Campus-Jun-2021_NS_20
## 1.3593685 1.9061547 0.9502705
## Campus-Jun-2021_S_1 Campus-Jun-2021_S_3 Campus-Jun-2021_S_5
## 2.2718685 2.0028830 1.7351265
## Campus-Jun-2021_S_7 Campus-Jun-2021_S_9 Campus-Jun-2021_S_11
## 0.9502705 2.1383972 1.5595812
## Campus-Jun-2021_S_13 Campus-Jun-2021_S_15 Campus-Jun-2021_S_17
## 1.8576598 0.7356219 2.1639557
## Campus-Jun-2021_S_19 Campus-maig-2020_NS_1 Campus-maig-2020_NS_2
## 1.9722470 0.8392696 1.3642812
## Campus-maig-2020_NS_3 Campus-maig-2020_NS_4 Campus-maig-2020_NS_5
## 1.6957425 2.5521719 1.5171064
## Campus-maig-2020_NS_6 Campus-maig-2020_NS_7 Campus-maig-2020_NS_8
## 2.2614951 1.2424533 1.7480673
## Campus-maig-2020_NS_9 Campus-maig-2020_NS_10 Campus-maig-2020_NS_11
## 1.6769878 1.9722470 1.7480673
## Campus-maig-2020_NS_12 Campus-maig-2020_NS_13 Campus-maig-2020_NS_14
## 1.0986123 1.5481155 1.8891592
## Campus-maig-2020_NS_15 Campus-maig-2020_NS_16 Campus-maig-2020_NS_17
## 1.8392967 1.0356097 1.3296613
## Campus-maig-2020_NS_18 Campus-maig-2020_NS_19 Campus-maig-2020_NS_20
## 1.5810938 1.3535914 0.9002561
## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 118 184.0553 24.55884 170.3457 9.341038 201.8116 140.7347 5.098803 81Corba i extrapolació per Franqueses:
#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades per Franqueses
abundacia_trampa_franqueses <- dades_finals %>%
filter(Zona == 'Franqueses') %>%
group_by(Mostreig, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxfranq <- cast(abundacia_trampa_franqueses, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxfranq <- as.data.frame(taula_auxfranq)
taula_auxfranq[is.na(taula_auxfranq)] <- 0
# Agrupar Mostreig amb Tractament i Trampa
data_shannon_franqueses <- taula_auxfranq %>%
mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament, -Trampa)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_franqueses) <- data_shannon_franqueses$Agrupacio
data_shannon_franqueses <- data_shannon_franqueses %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivfranq <- diversity(data_shannon_franqueses)
print(shannondivfranq)## Franqueses-Maig-2021_NS_1 Franqueses-Maig-2021_NS_2 Franqueses-Maig-2021_NS_3
## 1.0789922 0.8675632 0.6931472
## Franqueses-Maig-2021_NS_4 Franqueses-Maig-2021_NS_5 Franqueses-Maig-2021_S_1
## 2.2739657 0.7549968 1.3321790
## Franqueses-Maig-2021_S_2 Franqueses-Maig-2021_S_4
## 1.0986123 1.5607104
## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 26 52.32292 16.96108 42.625 8.605049 53.01786 33.13315 4.330973 8Corba i extrapolació per Sabadell
#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades per Sabadell
abundacia_trampa_sabadell <- dades_finals %>%
filter(Zona == 'Sabadell') %>%
group_by(Mostreig, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxsaba <- cast(abundacia_trampa_sabadell, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxsaba <- as.data.frame(taula_auxsaba)
taula_auxsaba[is.na(taula_auxsaba)] <- 0
# Agrupar Mostreig amb Tractament i Trampa
data_shannon_sabadell <- taula_auxsaba %>%
mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament, -Trampa)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_sabadell) <- data_shannon_sabadell$Agrupacio
data_shannon_sabadell <- data_shannon_sabadell %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivsaba <- diversity(data_shannon_sabadell)
print(shannondivsaba)## Sabadell-Juny-2021_NS_1 Sabadell-Juny-2021_NS_4 Sabadell-Juny-2021_NS_7
## 1.7478681 0.9502705 1.7917595
## Sabadell-Juny-2021_NS_8 Sabadell-Juny-2021_NS_11 Sabadell-Juny-2021_S_2
## 1.4750763 1.3862944 0.6365142
## Sabadell-Juny-2021_S_3 Sabadell-Juny-2021_S_5 Sabadell-Juny-2021_S_6
## 0.5623351 1.4995094 1.0986123
## Sabadell-Juny-2021_S_9 Sabadell-Juny-2021_S_10 Sabadell-Juny-2021_S_12
## 0.6931472 0.5004024 1.2852930
## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 23 45.45833 17.10731 35.83333 5.316379 44.4697 28.44918 2.631169 12T-student per comparar la diversitat de Shannon entre mostrejos i tractaments
#Convertir la taula de shannondiv en dataframe
shannon_ttest <- as.data.frame(shannondiv)
# Convertir l'ID de les files (mostreig + tractament + número de trampa) en una variable (nova columna)
shannon_ttest$row <- rownames(shannon_ttest)
# De la nova variable, separar en dues columnes el mostreig i eltractament (i no utilitzem el número de trampa)
shannon_ttest <- shannon_ttest %>% mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2], Mostreig = str_split(row, "_", simplify = TRUE)[,1]) %>% dplyr::select(-row)
# Prova t-student i boxplots
shannonT <- shannon_ttest %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(shannondiv[Tractament=='S']),
mean_NS = mean(shannondiv[Tractament=='NS']),
t_test_result = ifelse(length(unique(Tractament)) == 2,
t.test(shannondiv ~ Tractament)$p.value,
NA),.groups = "drop")
library(openxlsx)
write.xlsx(shannonT,"shannonT.xlsx")
ggplot(shannon_ttest, aes(x = factor(Mostreig), y = shannondiv, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Diversitat de Shannon x Zona i Tractaments",
x = "Mostreig",
y = "Diversitat Shannon") +
scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
scale_color_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
theme_minimal()
Només en el mostreig de Sabadell la diversitat de Shannon és significativament diferent (p valor < 0.05): l’índex de Shannon és major al tractament NS.
Diversitat de Shannon d’abelles (per mostreig)
# Preparar les dades: crear taula de dades d'abundància només per abelles
abundanciabelles_ID <- dades_finals %>%
mutate(Mostreig = paste(Zona, Mes, Any, sep = "-")) %>%
filter(Ordre == 'HymenopteraAbella') %>%
group_by(Mostreig, Tractament, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.# View(abundanciabelles_ID)
#Crear una taula auxiliar de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux <- cast(abundanciabelles_ID, Mostreig + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux <- as.data.frame(taula_aux)
taula_aux[is.na(taula_aux)] <- 0
#Agrupar Mostreig i tractament en una nova variable (Agrupacio)
data_shannon_abelles <- taula_aux %>%
mutate(Agrupacio = paste(Mostreig, Tractament, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament)
#Definir els rownames amb la variable Agrupacio (mostreig i tractament)
rownames(data_shannon_abelles) <- data_shannon_abelles$Agrupacio
data_shannon_abelles <- data_shannon_abelles %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
#View(data_shannon_abelles)# Càlcul diversitat de Shannon
shannondiv_abelles <- diversity(data_shannon_abelles)
print(shannondiv_abelles)## Campus-juliol-2020_NS Campus-juliol-2020_S Campus-Juliol-2021_NS
## 2.611229 1.905899 2.090031
## Campus-Juliol-2021_S Campus-Jun-2021_NS Campus-Jun-2021_S
## 1.806531 2.701045 2.080060
## Campus-maig-2020_NS Franqueses-Maig-2021_NS Franqueses-Maig-2021_S
## 2.199104 2.239746 1.560710
## Sabadell-Juny-2021_NS Sabadell-Juny-2021_S
## 2.045357 1.884985Shannon abelles per trampa
#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa
abelles_abundacia_trampa <- dades_finals %>%
filter(Ordre == 'HymenopteraAbella') %>%
group_by(Mostreig, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.## # A tibble: 371 × 5
## # Groups: Mostreig, Tractament, Trampa [100]
## Mostreig Tractament Trampa ID Numero_total
## <chr> <chr> <dbl> <chr> <int>
## 1 Campus-Juliol-2021 NS 2 Lasioglossum glabriusculum 1
## 2 Campus-Juliol-2021 NS 4 Lasioglossum malachurum 1
## 3 Campus-Juliol-2021 NS 6 Lasioglossum glabriusculum 1
## 4 Campus-Juliol-2021 NS 6 Lasioglossum politum 3
## 5 Campus-Juliol-2021 NS 8 Lasioglossum glabriusculum 3
## 6 Campus-Juliol-2021 NS 8 Megachile argentata 3
## 7 Campus-Juliol-2021 NS 8 Seladonia gemmea 1
## 8 Campus-Juliol-2021 NS 10 Bombus pascuorum 1
## 9 Campus-Juliol-2021 NS 12 Halictus sp. 1
## 10 Campus-Juliol-2021 NS 12 Lasioglossum glabriusculum 2
## # ℹ 361 more rows#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux_abelles <- cast(abelles_abundacia_trampa, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_aux_abelles <- as.data.frame(taula_aux_abelles)
taula_aux_abelles[is.na(taula_aux_abelles)] <- 0# Agrupar Mostreig amb Tractament i Trampa
abelles_shannon_MTTR <- taula_aux_abelles %>%
mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament, -Trampa)# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(abelles_shannon_MTTR) <- abelles_shannon_MTTR$Agrupacio
abelles_shannon_MTTR <- abelles_shannon_MTTR %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
abelles_shannondiv <- diversity(abelles_shannon_MTTR)
print(abelles_shannondiv)## Campus-juliol-2020_NS_2 Campus-juliol-2020_NS_3 Campus-juliol-2020_NS_4
## 1.3862944 1.6674619 1.6417347
## Campus-juliol-2020_NS_6 Campus-juliol-2020_NS_8 Campus-juliol-2020_NS_10
## 1.3296613 1.6674619 0.7963116
## Campus-juliol-2020_NS_11 Campus-juliol-2020_NS_12 Campus-juliol-2020_NS_13
## 0.0000000 1.5607104 1.6434177
## Campus-juliol-2020_NS_14 Campus-juliol-2020_NS_15 Campus-juliol-2020_NS_16
## 0.7549968 1.5498260 0.6931472
## Campus-juliol-2020_NS_20 Campus-juliol-2020_S_1 Campus-juliol-2020_S_5
## 1.7478681 0.6931472 1.3517840
## Campus-juliol-2020_S_7 Campus-juliol-2020_S_9 Campus-juliol-2020_S_17
## 0.5505734 0.0000000 1.0986123
## Campus-juliol-2020_S_18 Campus-juliol-2020_S_19 Campus-juliol-2020_S_21
## 1.3862944 0.6931472 1.1189689
## Campus-juliol-2020_S_22 Campus-Juliol-2021_NS_2 Campus-Juliol-2021_NS_4
## 0.6108643 0.0000000 0.0000000
## Campus-Juliol-2021_NS_6 Campus-Juliol-2021_NS_8 Campus-Juliol-2021_NS_10
## 0.5623351 1.0042425 0.0000000
## Campus-Juliol-2021_NS_12 Campus-Juliol-2021_NS_14 Campus-Juliol-2021_NS_16
## 1.0549202 1.0397208 0.0000000
## Campus-Juliol-2021_NS_18 Campus-Juliol-2021_NS_20 Campus-Juliol-2021_S_1
## 1.3862944 0.6931472 0.6931472
## Campus-Juliol-2021_S_3 Campus-Juliol-2021_S_5 Campus-Juliol-2021_S_7
## 0.6365142 1.8662160 1.4369665
## Campus-Juliol-2021_S_9 Campus-Juliol-2021_S_11 Campus-Juliol-2021_S_13
## 1.1973401 1.4388830 0.0000000
## Campus-Juliol-2021_S_15 Campus-Juliol-2021_S_17 Campus-Jun-2021_NS_2
## 0.6931472 0.5091373 1.4735024
## Campus-Jun-2021_NS_4 Campus-Jun-2021_NS_6 Campus-Jun-2021_NS_8
## 1.0986123 0.8675632 2.3393717
## Campus-Jun-2021_NS_10 Campus-Jun-2021_NS_12 Campus-Jun-2021_NS_14
## 1.3862944 2.0947290 1.4925477
## Campus-Jun-2021_NS_16 Campus-Jun-2021_NS_18 Campus-Jun-2021_NS_20
## 1.2130076 1.5607104 0.0000000
## Campus-Jun-2021_S_1 Campus-Jun-2021_S_3 Campus-Jun-2021_S_5
## 1.6094379 1.5498260 1.0397208
## Campus-Jun-2021_S_7 Campus-Jun-2021_S_9 Campus-Jun-2021_S_11
## 0.5623351 1.5810938 1.5595812
## Campus-Jun-2021_S_13 Campus-Jun-2021_S_15 Campus-Jun-2021_S_17
## 1.4388830 0.4101163 1.3321790
## Campus-Jun-2021_S_19 Campus-maig-2020_NS_1 Campus-maig-2020_NS_2
## 1.2424533 0.6682485 0.9611277
## Campus-maig-2020_NS_3 Campus-maig-2020_NS_4 Campus-maig-2020_NS_5
## 0.6365142 1.6094379 0.6365142
## Campus-maig-2020_NS_6 Campus-maig-2020_NS_7 Campus-maig-2020_NS_8
## 1.8845210 0.6931472 1.3862944
## Campus-maig-2020_NS_9 Campus-maig-2020_NS_10 Campus-maig-2020_NS_11
## 1.3321790 0.9502705 1.3862944
## Campus-maig-2020_NS_12 Campus-maig-2020_NS_13 Campus-maig-2020_NS_14
## 0.6931472 1.3114313 1.3296613
## Campus-maig-2020_NS_15 Campus-maig-2020_NS_16 Campus-maig-2020_NS_17
## 1.4978661 0.9502705 0.6365142
## Campus-maig-2020_NS_18 Campus-maig-2020_NS_19 Franqueses-Maig-2021_NS_1
## 1.3862944 0.9251291 0.6931472
## Franqueses-Maig-2021_NS_2 Franqueses-Maig-2021_NS_3 Franqueses-Maig-2021_NS_4
## 0.5004024 0.6931472 1.6674619
## Franqueses-Maig-2021_NS_5 Franqueses-Maig-2021_S_1 Franqueses-Maig-2021_S_2
## 1.0986123 0.0000000 1.0986123
## Franqueses-Maig-2021_S_4 Sabadell-Juny-2021_NS_1 Sabadell-Juny-2021_NS_4
## 0.6931472 1.7478681 0.0000000
## Sabadell-Juny-2021_NS_7 Sabadell-Juny-2021_NS_8 Sabadell-Juny-2021_NS_11
## 1.3862944 1.2424533 1.3862944
## Sabadell-Juny-2021_S_2 Sabadell-Juny-2021_S_3 Sabadell-Juny-2021_S_5
## 0.6365142 0.5623351 1.4995094
## Sabadell-Juny-2021_S_6 Sabadell-Juny-2021_S_9 Sabadell-Juny-2021_S_10
## 0.6931472 0.6931472 0.5004024
## Sabadell-Juny-2021_S_12
## 1.0986123#Convertir la taula de shannondiv en dataframe
abelles_shannon_ttest <- as.data.frame(abelles_shannondiv)
# Convertir l'ID de les files (mostreig + tractament + número de trampa) en una variable (nova columna)
abelles_shannon_ttest$row <- rownames(abelles_shannon_ttest)
# De la nova variable, separar en dues columnes el mostreig i eltractament (i no utilitzem el número de trampa)
abelles_shannon_ttest <- abelles_shannon_ttest %>% mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2], Mostreig = str_split(row, "_", simplify = TRUE)[,1]) %>% dplyr::select(-row)
# Prova t-student i boxplots
abelles_shannonT <- abelles_shannon_ttest %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(abelles_shannondiv[Tractament=='S']),
mean_NS = mean(abelles_shannondiv[Tractament=='NS']),
t_test_result = ifelse(length(unique(Tractament)) == 2,
t.test(abelles_shannondiv ~ Tractament)$p.value,
NA),.groups = "drop")
write.xlsx(abelles_shannonT,"abelles_shannonT.xlsx")
ggplot(abelles_shannon_ttest, aes(x = factor(Mostreig), y = abelles_shannondiv, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Diversitat de Shannon d'abelles x Zona i Tractaments",
x = "Mostreig",
y = "Diversitat Shannon") +
scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
scale_color_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
theme_minimal()
Corbes d’acumulació d’espècies per mostreig
Corba d’acumulació. Riquesa en funció del nº de trampes (specaccum)
# A la taula data_shannon_MTTR extreure una nova columna que indiqui el mostreig
metadata_shannon <- data_shannon_MTTR %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1])
# View(metadata_shannon)
# Subset each mostreig into its own dataframe (i eliminar la columna de "mostreig" perquè totes les dades de la taula siguin numèriques)
metadata_shannon %>% filter(Zona == 'Campus') %>% dplyr::select(-Zona)-> Campus
metadata_shannon %>% filter(Zona == 'Franqueses') %>% dplyr::select(-Zona) -> Franqueses
metadata_shannon %>% filter(Zona == 'Llig') %>% dplyr::select(-Zona) -> Llig
metadata_shannon %>% filter(Zona == 'Mura') %>% dplyr::select(-Zona) -> Mura
metadata_shannon %>% filter(Zona == 'Sabadell') %>% dplyr::select(-Zona) -> Sabadell
metadata_shannon %>% filter(Zona == 'TV') %>% dplyr::select(-Zona) -> TV7041
# Corba d'acumulació d'espècies per cada mostreig
corba1 = specaccum(Campus, method = "random")
corba2 = specaccum(Franqueses, method = "random")
corba3 = specaccum(Sabadell, method = "random")
# Creating data frames for each set of data
data1 <- data.frame(Sites = corba1$sites, Richness = corba1$richness, SD = corba1$sd)
data1$Zona <- "Campus"
data2 <- data.frame(Sites = corba2$sites, Richness = corba2$richness, SD = corba2$sd)
data2$Zona <- "Franqueses"
data3 <- data.frame(Sites = corba3$sites, Richness = corba3$richness, SD = corba3$sd)
data3$Zona <- "Sabadell"
# Combine all data frames into one
dades_combinades <- rbind(data1, data2, data3)
library(ggplot2)
# Gráfico de dispersión con líneas y error bars
ggplot(dades_combinades, aes(x = Sites, y = Richness, color = Zona)) +
geom_point() +
geom_line(aes(group = Zona), alpha = 0.7) + # Líneas
# geom_ribbon(aes(ymin = Richness - 2 * SD, ymax = Richness + 2 * SD, fill = Mostreig), alpha = 0.01) + # Incertidumbre
geom_errorbar(aes(ymin = Richness - 2 * SD, ymax = Richness + 2 * SD), width = 0.2, alpha = 0.5) + # Barras de error
ylim(0,NA)+
labs(x = "Nº trampes", y = "Riquesa", color = "Zona") +
theme_minimal()
Corba de rarefracció. Riquesa en funció del nombre d’individus (rarefy)
dades_agrupades <- metadata_shannon %>%
group_by(Zona) %>%
summarise(across(everything(), sum)) %>%
as.data.frame()
#View (dades_agrupades)
# Buscar el nombre mínim d'observacions en una zona (és 84, a TV juliol)
abundacia_trampa <- dades_finals %>%
group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.## # A tibble: 564 × 7
## # Groups: Zona, Mes, Any, Tractament, Trampa [101]
## Zona Mes Any Tractament Trampa ID Numero_total
## <chr> <chr> <dbl> <chr> <dbl> <chr> <int>
## 1 Campus Juliol 2021 NS 2 Lasioglossum glabriusculum 1
## 2 Campus Juliol 2021 NS 4 Curculionidae 1
## 3 Campus Juliol 2021 NS 4 Lasioglossum malachurum 1
## 4 Campus Juliol 2021 NS 4 Pompilidae 1
## 5 Campus Juliol 2021 NS 4 Pyronia cecilia 1
## 6 Campus Juliol 2021 NS 6 Lasioglossum glabriusculum 1
## 7 Campus Juliol 2021 NS 6 Lasioglossum politum 3
## 8 Campus Juliol 2021 NS 6 Pompilidae 1
## 9 Campus Juliol 2021 NS 8 Chrysididae 1
## 10 Campus Juliol 2021 NS 8 Hylaeus sp. 1
## # ℹ 554 more rowsmin_n <- abundacia_trampa %>%
group_by(Zona) %>%
summarize (n_especimens = sum(Numero_total)) %>%
summarize(min=min(n_especimens)) %>%
pull(min)
rownames(dades_agrupades) <- dades_agrupades$Zona
dades_agrupades <- dades_agrupades[,-1]
dades_agrupades <- dades_agrupades %>%
mutate_all(as.numeric)
vegan::rarefy(dades_agrupades,sample=min_n)## Campus Franqueses Sabadell
## 28.12144 26.00000 19.50913
## attr(,"Subsample")
## [1] 59col <- c("#81b29a","#3d405b", "#f2cc8f")
rarecurve(dades_agrupades, step = 5, col = col, lwd = 7, cex = 0.01, cex.axis = 2, cex.lab = 2)
legend("right", legend = rownames(dades_agrupades), col = col, lwd = 5, cex = 2, bty = "n")
Corba d’interpolació: riquesa en funció del nombre d’individus (iNEXT). (??)
library(iNEXT)
transposades <- t(dades_agrupades)
D_abund <- iNEXT (transposades, datatype = 'abundance')
plot (D_abund)
Corbes de rarefracció d’abelles per cada tractament
# Extreure la columna de tractament a la taula d'abundàncies per ID de cada mostreig
dades_shannon_abelles <- data_shannon_abelles %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2])
# Agrupar per tractament
dades_agrupades_tractament <- dades_shannon_abelles %>%
group_by(Tractament) %>%
summarise(across(everything(), sum)) %>%
as.data.frame()
#View (dades_agrupades_tractament)
dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
# Buscar el nombre mínim d'observacions en una zona
abundacia_trampa_abelles <- dades_abelles %>%
group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.#View (abundacia_trampa_abelles)
min_n <- abundacia_trampa_abelles %>%
group_by(Tractament) %>%
summarize (n_especimens = sum(Numero_total)) %>%
summarize(min=min(n_especimens)) %>%
pull(min)
rownames(dades_agrupades_tractament) <- dades_agrupades_tractament$Tractament
dades_agrupades_tractament <- dades_agrupades_tractament[,-1]
dades_agrupades_tractament <- dades_agrupades_tractament %>%
mutate_all(as.numeric)
vegan::rarefy(dades_agrupades_tractament,sample=min_n)## NS S
## 59.89088 37.00000
## attr(,"Subsample")
## [1] 321col <- c("#81b29a","#e07a5f")
rarecurve(dades_agrupades_tractament, step = 5, col = col, lwd = 2, cex = 0.01, cex.axis = 1.1, cex.lab = 1.1)
legend("right", legend = rownames(dades_agrupades_tractament), col = col, lwd = 3, cex = 1.1, bty = "n")
Corbes de rarefracció d’abelles al Campus per tractament
# Extreure la columna Zona de dades_shannon_abelles i filtrar per Campus
dades_shannon_abellescampus <- dades_shannon_abelles %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1])
dades_shannon_abellescampus <- dades_shannon_abellescampus[dades_shannon_abellescampus$Zona == "Campus", ]
# Agrupar per tractament
dadescampus_agrupades_tractament <- dades_shannon_abellescampus %>%
group_by(Tractament) %>%
summarise(across(where(is.numeric), sum)) %>%
as.data.frame()
#View (dadescampus_agrupades_tractament)
# Buscar el nombre mínim d'observacions en una zona
abundacia_trampa_abellescampus <- dades_abelles %>%
group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.#View (abundacia_trampa_abellescampus)
min_n <- abundacia_trampa_abellescampus %>%
group_by(Tractament) %>%
summarize (n_especimens = sum(Numero_total)) %>%
summarize(min=min(n_especimens)) %>%
pull(min)
rownames(dadescampus_agrupades_tractament) <- dadescampus_agrupades_tractament$Tractament
dadescampus_agrupades_tractament <- dadescampus_agrupades_tractament[,-1]
dadescampus_agrupades_tractament <- dadescampus_agrupades_tractament %>%
mutate_all(as.numeric)
vegan::rarefy(dadescampus_agrupades_tractament,sample=min_n)## Warning in vegan::rarefy(dadescampus_agrupades_tractament, sample = min_n):
## requested 'sample' was larger than smallest site maximum (266)## NS S
## 54.01949 35.00000
## attr(,"Subsample")
## [1] 321col <- c("#81b29a","#e07a5f")
rarecurve(dadescampus_agrupades_tractament, step = 5, col = col, lwd = 4, cex = 0.01, cex.axis = 1.1, cex.lab = 1.1)
legend("right", legend = rownames(dadescampus_agrupades_tractament), col = col, lwd = 3, cex = 1.1, bty = "n")
PROVES T-STUDENT
Comparació de la riquesa per mostrejos i tractaments
Riquesa total
riquesa_x_trampa <- dades_finals %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Riquesa = n_distinct(ID),.groups = "drop")
proves_t <- riquesa_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Riquesa[Tractament == "S"]),
mean_NS = mean(Riquesa[Tractament == "NS"]),
t_test_result = ifelse(length(unique(Tractament)) == 2,
t.test(Riquesa ~ Tractament)$p.value,
NA)
,.groups = "drop")
proves_t## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 4.33 3.1 0.303
## 2 Campus-Jun-2021 7.3 9.5 0.190
## 3 Campus-juliol-2020 3.89 5.92 0.0281
## 4 Campus-maig-2020 NaN 6.55 NA
## 5 Franqueses-Maig-2021 4 4.6 0.743
## 6 Sabadell-Juny-2021 3.43 4.8 0.201ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA, aes(color = Tractament)) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Riquesa per Zona i Tractaments",
x = "Mostreig",
y = "Riquesa") +
scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
scale_color_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
theme_minimal()
Generalment es compleix un patró de major riquesa al tractament NS (excepte a Campus-Juliol-2021). Al mostreig del Campus-Juliol_2020 aquesta diferència és significativa (p-valor= 0,0281)
Riquesa d’abelles
# Filtrar les dades d'abelles
dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
head(dades_abelles) ## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 159 Campus Juliol 7 2021 10 NS <NA>
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr># Proves t-student i plots
riquesa_x_trampa <- dades_abelles %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Riquesa = n_distinct(ID),.groups = "drop")
proves_t_abelles <- riquesa_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Riquesa[Tractament == "S"]),
mean_NS = mean(Riquesa[Tractament == "NS"]),
t_test_result = ifelse(length(unique(Tractament)) == 2,
t.test(Riquesa ~ Tractament)$p.value,
NA)
,.groups = "drop")
proves_t_abelles## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 3.89 2.1 0.0847
## 2 Campus-Jun-2021 4.3 5.4 0.338
## 3 Campus-juliol-2020 2.78 4.46 0.00961
## 4 Campus-maig-2020 NaN 3.84 NA
## 5 Franqueses-Maig-2021 2 3 0.341
## 6 Sabadell-Juny-2021 3.14 3.8 0.568ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Riquesa d'abelles per Zona i Tractaments",
x = "Mostreig",
y = "Riquesa") +
scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
scale_color_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
theme_minimal()
Generalment es compleix un patró de major riquesa d’abelles al tractament NS (excepte a Campus-Juliol-2021). Al mostreig del Campus-Juliol_2020 aquesta diferència és significativa (p-valor= 0.0096)
Riquesa de coleòpters
# Filtrar les dades per coleòpters
dades_coleoptera <- dades_finals[dades_finals$Ordre == "Coleoptera", ]
head(dades_coleoptera) ## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 29 Campus Juliol 7 2021 4 NS <NA>
## 2 158 Campus Juliol 7 2021 10 NS <NA>
## 3 183 Campus Juliol 7 2021 5 S <NA>
## 4 101 Franqueses Maig 5 2021 2 NS <NA>
## 5 150 Franqueses Maig 5 2021 5 NS Chiqui bupresti…
## 6 151 Franqueses Maig 5 2021 5 NS Chiqui bupresti…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr># Proves t-student i plots
riquesa_x_trampa <- dades_coleoptera %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Riquesa = n_distinct(ID),.groups = "drop")
proves_t_coleoptera <- riquesa_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Riquesa[Tractament == "S"]),
mean_NS = mean(Riquesa[Tractament == "NS"]),
t_test_result = tryCatch(
t.test(Riquesa ~ Tractament)$p.value,
error = function(e) NA
)
,.groups = "drop")
proves_t_coleoptera## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 1 1 NA
## 2 Campus-Jun-2021 2.88 3.7 0.357
## 3 Campus-juliol-2020 1.33 1.25 0.851
## 4 Campus-maig-2020 NaN 2.44 NA
## 5 Franqueses-Maig-2021 1.5 1 0.500
## 6 Sabadell-Juny-2021 1 1.33 0.423 ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Riquesa x Zona i Tractaments",
x = "Mostreig",
y = "Riquesa")
ERROR: NOT ENOUGH OBSERVATIONS. ALGUNS MOSTREJOS NO TENEN COLEÒPTERS.
Riquesa de sírfids
# Filtrar les dades per sírfids
dades_sirfids <-dades_finals %>% filter(Familia=="Syrphidae")
head(dades_sirfids) ## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 226 Campus Juliol 7 2021 7 S <NA>
## 2 265 Campus Juliol 7 2021 14 NS <NA>
## 3 277 Franqueses Maig 5 2021 4 NS <NA>
## 4 278 Franqueses Maig 5 2021 4 NS <NA>
## 5 279 Franqueses Maig 5 2021 4 NS <NA>
## 6 280 Franqueses Maig 5 2021 4 NS Episyrphus balt…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr># Proves t-student i plots
riquesa_x_trampa <- dades_sirfids %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Riquesa = n_distinct(ID),.groups = "drop")
proves_t_sirfids <- riquesa_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Riquesa[Tractament == "S"]),
mean_NS = mean(Riquesa[Tractament == "NS"]),
t_test_result = tryCatch(
t.test(Riquesa ~ Tractament)$p.value,
error = function(e) NA
)
,.groups = "drop")
proves_t_sirfids## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <lgl>
## 1 Campus-Juliol-2021 1 1 NA
## 2 Campus-Jun-2021 1 1 NA
## 3 Campus-juliol-2020 3 1 NA
## 4 Campus-maig-2020 NaN 1.18 NA
## 5 Franqueses-Maig-2021 1 3 NA
## 6 Sabadell-Juny-2021 NaN 1 NAggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Riquesa x Zona i Tractaments",
x = "Mostreig",
y = "Riquesa")
Hi ha poques dades per obtenir resultats del t-test.
Comparació de l’abundància per mostrejos i tractaments
Abundància total
abundancia_x_trampa <- dades_finals %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Abundancia = n(),.groups = "drop")
proves_t_abtotal <- abundancia_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Abundancia[Tractament == "S"]),
mean_NS = mean(Abundancia[Tractament == "NS"]),
t_test_result = ifelse(length(unique(Tractament)) == 2,
t.test(Abundancia ~ Tractament)$p.value,
NA)
,.groups = "drop")
proves_t_abtotal## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 14.1 4.5 0.130
## 2 Campus-Jun-2021 11.6 18.5 0.113
## 3 Campus-juliol-2020 8.89 8.54 0.904
## 4 Campus-maig-2020 NaN 15.4 NA
## 5 Franqueses-Maig-2021 4.67 9 0.178
## 6 Sabadell-Juny-2021 7.29 5.8 0.606 ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Abundancia x Zona i Tractaments",
x = "Mostreig",
y = "Abundancia") +
scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
scale_color_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
theme_minimal()
No hi ha cap resultat significatiu (p-valor<0.05). No hi ha un patró clar de com varia l’abundància en funció del tractament.
Abundància d’abelles
# Filtrar les dades d'abelles
dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
head(dades_abelles) ## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 159 Campus Juliol 7 2021 10 NS <NA>
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr># Proves t-student i plots
abundancia_x_trampa <- dades_abelles %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Abundancia = n(),.groups = "drop")
proves_t_ababelles <- abundancia_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Abundancia[Tractament == "S"]),
mean_NS = mean(Abundancia[Tractament == "NS"]),
t_test_result = ifelse(length(unique(Tractament)) == 2,
t.test(Abundancia ~ Tractament)$p.value,
NA)
,.groups = "drop")
proves_t_ababelles## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 13.7 3.5 0.108
## 2 Campus-Jun-2021 7.6 8.9 0.589
## 3 Campus-juliol-2020 7.44 7.08 0.896
## 4 Campus-maig-2020 NaN 8.68 NA
## 5 Franqueses-Maig-2021 2 4.4 0.0902
## 6 Sabadell-Juny-2021 7 4.8 0.454 ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Abundancia x Zona i Tractaments",
x = "Mostreig",
y = "Abundancia") +
scale_fill_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
scale_color_manual(values = c("S" = "#e07a5f", "NS" = "#81b29a")) +
theme_minimal()
A Franqueses l’abundància d’abelles és significativament superior al tractament NS (p-valor=0.09015). En la resta de mostrejos el patró no és clar.
Abundància de coleòpters
# Filtrar les dades per coleòpters
dades_coleoptera <- dades_finals[dades_finals$Ordre == "Coleoptera", ]
head(dades_coleoptera) ## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 29 Campus Juliol 7 2021 4 NS <NA>
## 2 158 Campus Juliol 7 2021 10 NS <NA>
## 3 183 Campus Juliol 7 2021 5 S <NA>
## 4 101 Franqueses Maig 5 2021 2 NS <NA>
## 5 150 Franqueses Maig 5 2021 5 NS Chiqui bupresti…
## 6 151 Franqueses Maig 5 2021 5 NS Chiqui bupresti…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr># Proves t-student i plots
abundancia_x_trampa <- dades_coleoptera %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Abundancia = n(),.groups = "drop")
proves_t_abcoleoptera <- abundancia_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Abundancia[Tractament == "S"]),
mean_NS = mean(Abundancia[Tractament == "NS"]),
t_test_result = tryCatch(
t.test(Abundancia ~ Tractament)$p.value,
error = function(e) NA
)
,.groups = "drop")
proves_t_abcoleoptera## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 1 1 NA
## 2 Campus-Jun-2021 4.12 9.2 0.175
## 3 Campus-juliol-2020 1.33 1.25 0.851
## 4 Campus-maig-2020 NaN 6.61 NA
## 5 Franqueses-Maig-2021 2 4.67 0.496
## 6 Sabadell-Juny-2021 1 1.33 0.423 ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Abundancia x Zona i Tractaments",
x = "Mostreig",
y = "Abundancia")
ERROR: NOT ENOUGH OBSERVATIONS. ALGUNS MOSTREJOS NO TENEN
COLEÒPTERS.
Abundància de sírfids
# Filtrar les dades per sírfids
dades_sirfids <-dades_finals %>% filter(Familia=="Syrphidae")
head(dades_sirfids) ## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 226 Campus Juliol 7 2021 7 S <NA>
## 2 265 Campus Juliol 7 2021 14 NS <NA>
## 3 277 Franqueses Maig 5 2021 4 NS <NA>
## 4 278 Franqueses Maig 5 2021 4 NS <NA>
## 5 279 Franqueses Maig 5 2021 4 NS <NA>
## 6 280 Franqueses Maig 5 2021 4 NS Episyrphus balt…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr># Proves t-student i plots
abundancia_x_trampa <- dades_sirfids %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Abundancia = n_distinct(ID),.groups = "drop")
proves_t_absirfids <- abundancia_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Abundancia[Tractament == "S"]),
mean_NS = mean(Abundancia[Tractament == "NS"]),
t_test_result = tryCatch(
t.test(Abundancia ~ Tractament)$p.value,
error = function(e) NA
)
,.groups = "drop")
proves_t_absirfids## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <lgl>
## 1 Campus-Juliol-2021 1 1 NA
## 2 Campus-Jun-2021 1 1 NA
## 3 Campus-juliol-2020 3 1 NA
## 4 Campus-maig-2020 NaN 1.18 NA
## 5 Franqueses-Maig-2021 1 3 NA
## 6 Sabadell-Juny-2021 NaN 1 NA ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Abundancia x Zona i Tractaments",
x = "Mostreig",
y = "Abundancia")
No hi ha prous dades per obtenir bons resultats de la prova t-student.
DIAGRAMES D’EULER:
Diagrames d’Euler totals per cada zona
# Preparació de les dades: càlcul de l'abundància de cada espècie (ID) a cada mostreig i tractament (NS, S i intersecció)
taula_aux_euler <- as.data.frame.matrix(table(interaction(dades_finals$Zona,dades_finals$ID),dades_finals$Tractament)) %>% mutate(`NS&S`=pmin(NS,S))
head(taula_aux_euler)## NS S NS&S
## Campus.Acmaeodera cylindrica 5 1 1
## Franqueses.Acmaeodera cylindrica 0 0 0
## Sabadell.Acmaeodera cylindrica 0 0 0
## Campus.Amegilla sp. 1 0 0
## Franqueses.Amegilla sp. 0 0 0
## Sabadell.Amegilla sp. 0 0 0# Reorganitzar la taula auxiliar perquè mostri el nom de les espècies (ID) trobades en cada mostreig i tractament
taula_euler_especies <- taula_aux_euler %>% mutate(Especie = sub("^[^.]+\\.", "", rownames(.)), Zona = sub("\\..*","",rownames(taula_aux_euler)))
taula_euler_especies <- taula_euler_especies %>% group_by(Zona) %>% summarise( S = list(unique(Especie[ S > 0 ])),
NS = list(unique(Especie[ NS > 0 ])),
`NS&S` = list(unique(Especie[ `NS&S` > 0 ])) )
head(taula_euler_especies)## # A tibble: 3 × 4
## Zona S NS `NS&S`
## <chr> <list> <list> <list>
## 1 Campus <chr [59]> <chr [98]> <chr [39]>
## 2 Franqueses <chr [11]> <chr [20]> <chr [5]>
## 3 Sabadell <chr [12]> <chr [16]> <chr [5]># Descarregar les llistes d'espècies en un excel
write.xlsx(taula_euler_especies, file.path(directorio, "2.0_Euler_Especies.xlsx"))# Agrupació per mostrejos (Zona, Mes, Any)
euler_data <- taula_euler_especies %>%
mutate(NS = map_dbl(NS, length)-map_dbl(`NS&S`, length),
S = map_dbl(S, length)-map_dbl(`NS&S`, length),
`NS&S` = map_dbl(`NS&S`, length))
head(euler_data)## # A tibble: 3 × 4
## Zona S NS `NS&S`
## <chr> <dbl> <dbl> <dbl>
## 1 Campus 20 59 39
## 2 Franqueses 6 15 5
## 3 Sabadell 7 11 5# Creació del diagrama d'Euler
apply_euler <- function(row){
euler(c("NS"=row$NS,"S"=row$S,"NS&S"=row$`NS&S`))
}
euler_data <- euler_data %>% rowwise() %>% mutate(euler_result= list(apply_euler(cur_data())))# Plot del diagrama d'Euler
quantities_font <- list(color = "white", size = 1.8, style = 2)
generate_euler_plots <- function(euler_result, zona, colors, label_size = 1.5){
plot(euler_result, main = paste(zona), quantities = list(type = c("counts"),
col = quantities_font$color,
cex = quantities_font$size,
font = quantities_font$style), fill = colors,
labels = list(col = "white", cex = label_size))
}
# Example of specifying a color palette
colors_palette <- c("#81b29a", "#e07a5f", "#3D405B")
# Use pmap to apply function with color argument
pmap(list(euler_data$euler_result, euler_data$Zona, list(colors_palette)), generate_euler_plots)## [[1]]
##
## [[2]]
##
## [[3]]
Diagrames d’Euler d’abelles (amb dades_abelles) per cada zona
dades_abelles <- dades_abelles[!(dades_abelles$Mes == "maig" & dades_abelles$Any == 2020 & dades_abelles$Zona == "Campus"), ]
# Preparació de les dades: càlcul de l'abundància de cada espècie (ID) a cada mostreig i tractament (NS, S i intersecció)
taula_aux_eulerabelles <- as.data.frame.matrix(table(interaction(dades_abelles$Zona,dades_abelles$ID),dades_abelles$Tractament)) %>% mutate(`NS&S`=pmin(NS,S))
head(taula_aux_eulerabelles)## NS S NS&S
## Campus.Amegilla sp. 1 0 0
## Franqueses.Amegilla sp. 0 0 0
## Sabadell.Amegilla sp. 0 0 0
## Campus.Andrena cinerea 0 0 0
## Franqueses.Andrena cinerea 1 0 0
## Sabadell.Andrena cinerea 0 0 0# Reorganitzar la taula auxiliar perquè mostri el nom de les espècies (ID) trobades en cada mostreig i tractament
taula_euler_abelles <- taula_aux_eulerabelles %>% mutate(Especie = sub("^[^.]+\\.", "", rownames(.)), Zona = sub("\\..*","",rownames(taula_aux_eulerabelles)))
taula_euler_abelles <- taula_euler_abelles %>% group_by(Zona) %>% summarise( S = list(unique(Especie[ S > 0 ])),
NS = list(unique(Especie[ NS > 0 ])),
`NS&S` = list(unique(Especie[ `NS&S` > 0 ])) )
head(taula_euler_abelles)## # A tibble: 3 × 4
## Zona S NS `NS&S`
## <chr> <list> <list> <list>
## 1 Campus <chr [35]> <chr [47]> <chr [21]>
## 2 Franqueses <chr [5]> <chr [12]> <chr [3]>
## 3 Sabadell <chr [10]> <chr [12]> <chr [5]># Descarregar les llistes d'espècies en un excel
write.xlsx(taula_euler_abelles, file.path(directorio, "Euler_Abelles.xlsx"))# Agrupació per mostrejos (Zona, Mes, Any)
euler_data_abelles <- taula_euler_abelles %>%
mutate(NS = map_dbl(NS, length)-map_dbl(`NS&S`, length),
S = map_dbl(S, length)-map_dbl(`NS&S`, length),
`NS&S` = map_dbl(`NS&S`, length))
head(euler_data_abelles)## # A tibble: 3 × 4
## Zona S NS `NS&S`
## <chr> <dbl> <dbl> <dbl>
## 1 Campus 14 26 21
## 2 Franqueses 2 9 3
## 3 Sabadell 5 7 5# Creació del diagrama d'Euler
apply_euler <- function(row){
euler(c("NS"=row$NS,"S"=row$S,"NS&S"=row$`NS&S`))
}
euler_data_abelles <- euler_data_abelles %>% rowwise() %>% mutate(euler_result= list(apply_euler(cur_data())))# Plot del diagrama d'Euler
quantities_font <- list(color = "white", size = 1.8, style = 2)
generate_euler_plots <- function(euler_result, zona, colors, label_size = 1.5){
plot(euler_result, main = paste(zona), quantities = list(type = c("counts"),
col = quantities_font$color,
cex = quantities_font$size,
font = quantities_font$style), fill = colors,
labels = list(col = "white", cex = label_size))
}
# Example of specifying a color palette
colors_palette <- c("#81b29a", "#e07a5f", "#3D405B")
# Use pmap to apply function with color argument
pmap(list(euler_data_abelles$euler_result, euler_data$Zona, list(colors_palette)), generate_euler_plots)## [[1]]
##
## [[2]]
##
## [[3]]
ANÀLISIS DE COMPONENTS PRINCIPALS (PCA)
PCA total
# Preparar dades: taula de dades Shannon per mostreig i tractament (MT). Normalitzar les dades: nombre d'individus entre el nombre de pantraps del mostreig
pca_data <- data_shannon_MT
# Calcular el nombre de trampes per cada mostreig i tractament
n_trampes <- dades_finals %>% group_by(Mostreig, Tractament) %>% summarise(numtrampes=n_distinct(Trampa)) %>% ungroup() %>% mutate(row=paste(Mostreig,Tractament,sep="_")) %>% dplyr::select(-Mostreig,-Tractament)## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.# Unir les dues taules relacionant les files del mateix mostreig i tractament
pca_data$row <- rownames(pca_data)
pca_data <- pca_data %>% inner_join(n_trampes, by=join_by(row))
# Dividir les abundàncies de cada ID entre el nombre de trampes de cada mostreig i tractament
num_columns <-ncol(pca_data)
pca_data[, 1:(num_columns - 2)] <- pca_data[, 1:(num_columns - 2)] / pca_data[, num_columns]
# Eliminar l'última columna de la taula (numtrampes), ja no cal per fer la PCA
pca_data <- pca_data %>% dplyr::select(-numtrampes)
pca_data<- pca_data %>%
mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2])# Mostrar la PCA
ggbiplot(pca_comps, labels = pca_data$row, groups = pca_data$Tractament,
ellipse = TRUE, var.axes = FALSE) +
theme_minimal() +
theme(text = element_text(size = 20))
PCA abelles
# Preparar dades: taula de dades Shannon d'abelles. Normalitzar les dades: nombre d'individus entre el nombre de pantraps del mostreig
pca_data_abelles <- data_shannon_abelles
# Calcular el nombre de trampes per cada mostreig i tractament
n_trampes_abelles <- dades_finals %>%
filter(Ordre == 'HymenopteraAbella') %>%
group_by(Mostreig, Tractament) %>%
summarise(numtrampes = n_distinct(Trampa)) %>%
ungroup() %>%
mutate(row = paste(Mostreig, Tractament, sep = "_")) %>%
dplyr::select(-Mostreig, -Tractament)## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.# Unir les dues taules relacionant les files del mateix mostreig i tractament
pca_data_abelles$row <- rownames(pca_data_abelles)
pca_data_abelles <- pca_data_abelles %>% inner_join(n_trampes_abelles, by=join_by(row))
# Dividir les abundàncies de cada ID entre el nombre de trampes de cada mostreig i tractament
num_columnsabe <-ncol(pca_data_abelles)
pca_data_abelles[, 1:(num_columnsabe - 2)] <- pca_data_abelles[, 1:(num_columnsabe - 2)] / pca_data_abelles[, num_columnsabe]
# Eliminar l'última columna de la taula (numtrampes), ja no cal per fer la PCA
pca_data_abelles <- pca_data_abelles %>% dplyr::select(-numtrampes)
pca_data_abelles<- pca_data_abelles %>%
mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2])# Fer la PCA
pca_abelles <- prcomp(pca_data_abelles %>% dplyr::select(-row,-Tractament))
summary(pca_abelles)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.9981 1.1729 0.61490 0.52897 0.48323 0.38196 0.31303
## Proportion of Variance 0.6025 0.2076 0.05706 0.04223 0.03524 0.02202 0.01479
## Cumulative Proportion 0.6025 0.8101 0.86717 0.90940 0.94464 0.96665 0.98144
## PC8 PC9 PC10 PC11
## Standard deviation 0.2482 0.19318 0.15503 2.126e-16
## Proportion of Variance 0.0093 0.00563 0.00363 0.000e+00
## Cumulative Proportion 0.9907 0.99637 1.00000 1.000e+00# Mostrar la PCA
ggbiplot(pca_abelles, labels = pca_data_abelles$row, groups = pca_data_abelles$Tractament,
ellipse = TRUE, var.axes = FALSE) +
theme_minimal() +
theme(text = element_text(size = 20))
PCA abelles Campus (per trampes. No sé ben bé què he fet)
# Preparar dades: taula de dades Shannon d'abelles per mostreig, tractament i trampa (dades NO normalitzades)
pca_data_campus <- data_shannon_MTTR %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
filter(Zona == 'Campus') %>%
mutate(Tractament = str_extract(rownames(.), "(?<=_)[^_]+(?=_)"))
#View(pca_data_campus)
# Fer la PCA
pca_campus <- prcomp(pca_data_campus %>% dplyr::select(-Zona,-Tractament))
# Mostrar la PCA
ggbiplot(pca_campus, labels = pca_data_campus$row, groups = pca_data_campus$Tractament,
ellipse = TRUE, var.axes = FALSE) +
theme_minimal() +
theme(text = element_text(size = 20))
NMDS, ANÀLISIS DE SIMILITUD (ANOSIM), ADONIS, I ESPÈCIES INDICADORES
NMDS pol·linitzadors Campus
dataNMDS_polcampus <- pca_data_campus %>% dplyr::select(-Zona, -Tractament)
# Extreure de la taula el tractament de cada mostreig: una nova columna
Metadata1 <- dataNMDS_polcampus %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Metadata1 <- Metadata1 %>% mutate(color=ifelse(Tractament=="NS","black","black"))# NMDS Pol·linitzadors Campus:
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)
set.seed(7)
# k=2: nombre de dimensions
NMDS_polcampus <- metaMDS(dataNMDS_polcampus, k=2)## Wisconsin double standardization
## Run 0 stress 0.1731463
## Run 1 stress 0.1733912
## ... Procrustes: rmse 0.0970393 max resid 0.5241494
## Run 2 stress 0.1727394
## ... New best solution
## ... Procrustes: rmse 0.08088437 max resid 0.4655322
## Run 3 stress 0.1732186
## ... Procrustes: rmse 0.08060157 max resid 0.4245035
## Run 4 stress 0.173302
## Run 5 stress 0.1744142
## Run 6 stress 0.1742466
## Run 7 stress 0.1744684
## Run 8 stress 0.1740094
## Run 9 stress 0.1740286
## Run 10 stress 0.1740606
## Run 11 stress 0.1750397
## Run 12 stress 0.1743163
## Run 13 stress 0.1747583
## Run 14 stress 0.1748468
## Run 15 stress 0.1731277
## ... Procrustes: rmse 0.01610656 max resid 0.0950575
## Run 16 stress 0.1732682
## Run 17 stress 0.1754394
## Run 18 stress 0.1734799
## Run 19 stress 0.1737735
## Run 20 stress 0.1729431
## ... Procrustes: rmse 0.05032583 max resid 0.2448404
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 16: no. of iterations >= maxit
## 4: stress ratio > sratmaxcolors <- c("#74AA90", "#DB6443")
# Plot NMDS pol·linitzadors Campus
ordiplot(NMDS_polcampus,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata1$Tractament)
for (i in seq_along(unique_groups)) {
ordihull(NMDS_polcampus, groups = Metadata1$Tractament,
show.groups = unique_groups[i],
draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_polcampus,display="species",col="black",air=1, cex = 1.2)
# Save the last plot
g <- recordPlot()
# Replay and save as PNG
width <- 17 # Width in inches
height <- 14 # Height in inches
res <- 500 # Resolution in dpi
# Replay and save as PNG with higher resolution and larger size
png("PolCampus.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()## quartz_off_screen
## 2Anàlisi de similitud (ANOSIM) Total
##
## Call:
## anosim(x = dataNMDS_polcampus, grouping = Metadata1$Tractament, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.0432
## Significance: 0.1186
##
## Permutation: free
## Number of permutations: 9999No signiticatiu (p valor > 0.05): la composició de les comunitats no és significativament diferent entre els tractaments S i NS. També s’observa al plot de NMDS: els polígons que corresponen als dos tractaments estan superposats. Això probablement sigui per haver fet la prova amb totes les zones de mostreig alhora.
ADONIS total
# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_Polcampus <- dataNMDS_polcampus %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
NMDStot.dist <- vegdist(wisconsin(sqrt(dataNMDS_polcampus)), k=2)
adonis(formula = NMDStot.dist~Tractament, data = Tractament_metadata_Polcampus, permutations = 999)$aov.tab## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tractament 1 0.837 0.83702 1.9865 0.02453 0.004 **
## Residuals 79 33.286 0.42135 0.97547
## Total 80 34.123 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1No hi ha diferències significatives entre tractaments (p-valor = 0.25874), però la interacció entre zona i tractament explica un 45% de la variabilitat, i és un resultat significatiu (p-valor=0.03996).
Espècies indicadores dels tractaments S i NS total
indicadores = multipatt(dataNMDS_polcampus, Metadata1$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.05
##
## Total number of species: 137
## Selected number of species: 7
## Number of species associated to 1 group: 7
##
## List of species associated to each combination:
##
## Group NS #sps. 4
## stat p.value
## Lasioglossum cf. pauxillum 0.253 0.0356 *
## Rhodanthidium septemdentatum 0.249 0.0181 *
## Panurgus dentipes 0.235 0.0413 *
## Oedemera nobilis 0.217 0.0305 *
##
## Group S #sps. 3
## stat p.value
## Lasioglossum politum 0.339 0.0003 ***
## Lasioglossum glabriusculum 0.246 0.0124 *
## Lasioglossum griseolum 0.198 0.0408 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1NMDS abelles Campus
dataNMDS_abecampus <- abelles_shannon_MTTR %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
filter(Zona == "Campus") %>% dplyr::select(-Zona)
# Extreure de la taula el tractament de cada mostreig: una nova columna
Metadata2 <- dataNMDS_abecampus %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Metadata2 <- Metadata2 %>% mutate(color=ifelse(Tractament=="NS","black","black"))# NMDS Abelles Campus:
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)
set.seed(7)
# k=2: nombre de dimensions
NMDS_abecampus <- metaMDS(dataNMDS_abecampus, k=2)## Wisconsin double standardization
## Run 0 stress 0.1265607
## Run 1 stress 0.1259806
## ... New best solution
## ... Procrustes: rmse 0.04368297 max resid 0.2210331
## Run 2 stress 0.1252628
## ... New best solution
## ... Procrustes: rmse 0.05143489 max resid 0.2737022
## Run 3 stress 0.1253757
## ... Procrustes: rmse 0.03273285 max resid 0.1846694
## Run 4 stress 0.1245387
## ... New best solution
## ... Procrustes: rmse 0.04125358 max resid 0.2741621
## Run 5 stress 0.1252543
## Run 6 stress 0.1252991
## Run 7 stress 0.1251473
## Run 8 stress 0.1252565
## Run 9 stress 0.125563
## Run 10 stress 0.1258315
## Run 11 stress 0.1249854
## ... Procrustes: rmse 0.03920178 max resid 0.2852872
## Run 12 stress 0.125473
## Run 13 stress 0.1247876
## ... Procrustes: rmse 0.03278848 max resid 0.2426446
## Run 14 stress 0.1255484
## Run 15 stress 0.1256229
## Run 16 stress 0.1266907
## Run 17 stress 0.1244476
## ... New best solution
## ... Procrustes: rmse 0.02375557 max resid 0.1366641
## Run 18 stress 0.1251499
## Run 19 stress 0.125273
## Run 20 stress 0.1256626
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: no. of iterations >= maxitcolors <- c("#74AA90", "#DB6443")
# Plot NMDS Abelles Campus
ordiplot(NMDS_abecampus,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata2$Tractament)
for (i in seq_along(unique_groups)) {
ordihull(NMDS_abecampus, groups = Metadata2$Tractament,
show.groups = unique_groups[i],
draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_abecampus,display="species",col="black",air=1, cex = 1.2)
# Save the last plot
g <- recordPlot()
# Replay and save as PNG
width <- 17 # Width in inches
height <- 14 # Height in inches
res <- 500 # Resolution in dpi
# Replay and save as PNG with higher resolution and larger size
png("AbeCampus.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()## quartz_off_screen
## 2Anàlisi de similitud (ANOSIM) Abelles Campus
##
## Call:
## anosim(x = dataNMDS_abecampus, grouping = Metadata2$Tractament, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.029
## Significance: 0.1956
##
## Permutation: free
## Number of permutations: 9999No signiticatiu (p valor = 0.1956): la composició de les comunitats no és significativament diferent entre els tractaments S i NS.
ADONIS Abelles Campus
# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_Abecampus <- dataNMDS_abecampus %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
NMDStot.dist <- vegdist(wisconsin(sqrt(dataNMDS_abecampus)), k=2)
adonis(formula = NMDStot.dist~Tractament, data = Tractament_metadata_Abecampus, permutations = 999)$aov.tab## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tractament 1 0.880 0.87965 2.1685 0.02705 0.006 **
## Residuals 78 31.641 0.40566 0.97295
## Total 79 32.521 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El factor Tractament explica una pequeña proporción de la variabilidad total (R2 = 0.02705). La diferencia en la composición de las comunidades entre los grupos definidos por Tractament es estadísticamente significativa (valor p = 0.006). El estadístico F.Model de 2.1685 sugiere que hay una diferencia entre los grupos mayor que la esperada por el azar. Por lo tanto, a pesar de que la proporción de la varianza explicada es pequeña, la diferencia en la composición entre los grupos de Tractament es significativa desde un punto de vista estadístico.
Espècies indicadores abelles Campus
indicadores2 = multipatt(dataNMDS_abecampus, Metadata2$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores2)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.05
##
## Total number of species: 80
## Selected number of species: 6
## Number of species associated to 1 group: 6
##
## List of species associated to each combination:
##
## Group NS #sps. 3
## stat p.value
## Lasioglossum cf. pauxillum 0.256 0.0362 *
## Rhodanthidium septemdentatum 0.252 0.0146 *
## Panurgus dentipes 0.241 0.0377 *
##
## Group S #sps. 3
## stat p.value
## Lasioglossum politum 0.338 0.0005 ***
## Lasioglossum glabriusculum 0.244 0.0162 *
## Lasioglossum griseolum 0.197 0.0435 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sí que hi ha espècies indicadores!!!!
NMDS pol·linitzadors Franqueses
dataNMDS_polfranqueses <- data_shannon_MTTR %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
filter(Zona == "Franqueses") %>% dplyr::select(-Zona)
# Extreure de la taula el tractament de cada mostreig: una nova columna
Metadata3 <- dataNMDS_polfranqueses %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Metadata3 <- Metadata3 %>% mutate(color=ifelse(Tractament=="NS","black","black"))# NMDS:
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)
set.seed(7)
# k=2: nombre de dimensions
NMDS_polfranqueses <- metaMDS(dataNMDS_polfranqueses, k=2)## Wisconsin double standardization
## Run 0 stress 9.499321e-05
## Run 1 stress 1.701317e-05
## ... New best solution
## ... Procrustes: rmse 0.1971238 max resid 0.285099
## Run 2 stress 0.04788483
## Run 3 stress 8.301434e-05
## ... Procrustes: rmse 0.1625411 max resid 0.2817144
## Run 4 stress 9.344796e-05
## ... Procrustes: rmse 0.2418918 max resid 0.3678843
## Run 5 stress 8.989242e-05
## ... Procrustes: rmse 0.1362289 max resid 0.221904
## Run 6 stress 7.231338e-05
## ... Procrustes: rmse 0.2054776 max resid 0.4848038
## Run 7 stress 0
## ... New best solution
## ... Procrustes: rmse 0.1785951 max resid 0.3935142
## Run 8 stress 0
## ... Procrustes: rmse 0.2030769 max resid 0.2694547
## Run 9 stress 0.0477509
## Run 10 stress 9.984028e-05
## ... Procrustes: rmse 0.1689318 max resid 0.2881185
## Run 11 stress 9.886818e-05
## ... Procrustes: rmse 0.2345039 max resid 0.4475393
## Run 12 stress 8.805342e-05
## ... Procrustes: rmse 0.2059719 max resid 0.3315155
## Run 13 stress 9.00509e-05
## ... Procrustes: rmse 0.1985079 max resid 0.2915023
## Run 14 stress 0.06801037
## Run 15 stress 0.04833749
## Run 16 stress 8.45766e-05
## ... Procrustes: rmse 0.1123823 max resid 0.1615996
## Run 17 stress 0
## ... Procrustes: rmse 0.2343765 max resid 0.4118304
## Run 18 stress 0
## ... Procrustes: rmse 0.2541536 max resid 0.4202579
## Run 19 stress 9.431488e-05
## ... Procrustes: rmse 0.1992038 max resid 0.3022094
## Run 20 stress 0.04775453
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 4: no. of iterations >= maxit
## 15: stress < smin
## 1: stress ratio > sratmax## Warning in metaMDS(dataNMDS_polfranqueses, k = 2): stress is (nearly) zero: you
## may have insufficient data## Warning in postMDS(out$points, dis, plot = max(0, plot - 1), ...): skipping
## half-change scaling: too few points below thresholdcolors <- c("#74AA90", "#DB6443")
# Plot NMDS pol·linitzadors Campus
ordiplot(NMDS_polfranqueses,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata3$Tractament)
for (i in seq_along(unique_groups)) {
ordihull(NMDS_polfranqueses, groups = Metadata3$Tractament,
show.groups = unique_groups[i],
draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_polfranqueses,display="species",col="black",air=1, cex = 1.2)
# Save the last plot
g <- recordPlot()
# Replay and save as PNG
width <- 17 # Width in inches
height <- 14 # Height in inches
res <- 500 # Resolution in dpi
# Replay and save as PNG with higher resolution and larger size
png("PolFranqueses.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()## quartz_off_screen
## 2Anàlisi de similitud (ANOSIM) Pol·linitzadors Franqueses
ano = anosim(dataNMDS_polfranqueses, Metadata3$Tractament, distance = "bray", permutations = 9999)
ano##
## Call:
## anosim(x = dataNMDS_polfranqueses, grouping = Metadata3$Tractament, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: -0.005128
## Significance: 0.4773
##
## Permutation: free
## Number of permutations: 9999No signiticatiu (p valor = 0.4773): la composició de les comunitats no és significativament diferent entre els tractaments S i NS.
ADONIS Pol·linitzadors Franqueses
# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_polfranqueses <- dataNMDS_polfranqueses %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
head(Tractament_metadata_polfranqueses)## Tractament
## Franqueses-Maig-2021_NS_1 NS
## Franqueses-Maig-2021_NS_2 NS
## Franqueses-Maig-2021_NS_3 NS
## Franqueses-Maig-2021_NS_4 NS
## Franqueses-Maig-2021_NS_5 NS
## Franqueses-Maig-2021_S_1 SNMDStot.dist <- vegdist(wisconsin(sqrt(dataNMDS_polfranqueses)), k=2)
adonis(formula = NMDStot.dist~Tractament, data = Tractament_metadata_polfranqueses, permutations = 9999)$aov.tab## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").## Permutation: free
## Number of permutations: 9999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tractament 1 0.39194 0.39194 0.89501 0.12981 0.6102
## Residuals 6 2.62749 0.43791 0.87019
## Total 7 3.01943 1.00000No hi ha diferències significatives entre tractaments (p-valor = 0.6102).
Espècies indicadores Pol·linitzadors Franqueses
indicadores3 = multipatt(dataNMDS_polfranqueses, Metadata3$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores3)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.05
##
## Total number of species: 137
## Selected number of species: 0
## Number of species associated to 1 group: 0
##
## List of species associated to each combination:
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1No hi ha espècies indicadores a Franqueses
NMDS abelles Franqueses
dataNMDS_abefranqueses <- abelles_shannon_MTTR %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
filter(Zona == "Franqueses") %>% dplyr::select(-Zona)
# Extreure de la taula el tractament de cada mostreig: una nova columna
Metadata4 <- dataNMDS_abefranqueses %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Metadata4 <- Metadata4 %>% mutate(color=ifelse(Tractament=="NS","black","black"))# NMDS abelles Franqueses:
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)
set.seed(7)
# k=2: nombre de dimensions
NMDS_abefranqueses <- metaMDS(dataNMDS_abefranqueses, k=2)## Run 0 stress 0
## Run 1 stress 0
## ... Procrustes: rmse 0.1665426 max resid 0.251658
## Run 2 stress 0
## ... Procrustes: rmse 0.09615897 max resid 0.1451541
## Run 3 stress 0
## ... Procrustes: rmse 0.134433 max resid 0.2182935
## Run 4 stress 9.652733e-05
## ... Procrustes: rmse 0.1826024 max resid 0.2600087
## Run 5 stress 9.726088e-05
## ... Procrustes: rmse 0.1430138 max resid 0.2252209
## Run 6 stress 3.030201e-05
## ... Procrustes: rmse 0.1575005 max resid 0.2675314
## Run 7 stress 0
## ... Procrustes: rmse 0.1895275 max resid 0.3082301
## Run 8 stress 0
## ... Procrustes: rmse 0.117114 max resid 0.1693173
## Run 9 stress 0
## ... Procrustes: rmse 0.1365075 max resid 0.1737913
## Run 10 stress 0
## ... Procrustes: rmse 0.09004667 max resid 0.1559889
## Run 11 stress 0
## ... Procrustes: rmse 0.1075005 max resid 0.1820846
## Run 12 stress 5.509784e-05
## ... Procrustes: rmse 0.1815644 max resid 0.2208323
## Run 13 stress 9.48924e-05
## ... Procrustes: rmse 0.2171857 max resid 0.4526091
## Run 14 stress 5.958231e-05
## ... Procrustes: rmse 0.1243509 max resid 0.1690826
## Run 15 stress 9.850827e-05
## ... Procrustes: rmse 0.2171258 max resid 0.452336
## Run 16 stress 4.866928e-05
## ... Procrustes: rmse 0.2540161 max resid 0.507386
## Run 17 stress 0
## ... Procrustes: rmse 0.2351303 max resid 0.4343606
## Run 18 stress 0
## ... Procrustes: rmse 0.1837085 max resid 0.2811085
## Run 19 stress 0
## ... Procrustes: rmse 0.184819 max resid 0.3921269
## Run 20 stress 9.072194e-05
## ... Procrustes: rmse 0.1452401 max resid 0.1910837
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress < smin## Warning in metaMDS(dataNMDS_abefranqueses, k = 2): stress is (nearly) zero: you
## may have insufficient data## Warning in postMDS(out$points, dis, plot = max(0, plot - 1), ...): skipping
## half-change scaling: too few points below thresholdcolors <- c("#74AA90", "#DB6443")
# Plot NMDS pol·linitzadors Campus
ordiplot(NMDS_abefranqueses,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata4$Tractament)
for (i in seq_along(unique_groups)) {
ordihull(NMDS_abefranqueses, groups = Metadata4$Tractament,
show.groups = unique_groups[i],
draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_abefranqueses,display="species",col="black",air=1, cex = 1.2)
# Save the last plot
g <- recordPlot()
# Replay and save as PNG
width <- 17 # Width in inches
height <- 14 # Height in inches
res <- 500 # Resolution in dpi
# Replay and save as PNG with higher resolution and larger size
png("AbeFranqueses.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()## quartz_off_screen
## 2Anàlisi de similitud (ANOSIM) Abelles Franqueses
ano = anosim(dataNMDS_abefranqueses, Metadata4$Tractament, distance = "bray", permutations = 9999)
ano##
## Call:
## anosim(x = dataNMDS_abefranqueses, grouping = Metadata4$Tractament, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.1128
## Significance: 0.2478
##
## Permutation: free
## Number of permutations: 9999No signiticatiu (p valor = 0.2478): la composició de les comunitats no és significativament diferent entre els tractaments S i NS.
ADONIS Abelles Franqueses
# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_abefranqueses <- dataNMDS_abefranqueses %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
head(Tractament_metadata_abefranqueses)## Tractament
## Franqueses-Maig-2021_NS_1 NS
## Franqueses-Maig-2021_NS_2 NS
## Franqueses-Maig-2021_NS_3 NS
## Franqueses-Maig-2021_NS_4 NS
## Franqueses-Maig-2021_NS_5 NS
## Franqueses-Maig-2021_S_1 SNMDStot.dist <- vegdist(wisconsin(sqrt(dataNMDS_abefranqueses)), k=2)
adonis(formula = NMDStot.dist~Tractament, data = Tractament_metadata_abefranqueses, permutations = 9999)$aov.tab## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").## Permutation: free
## Number of permutations: 9999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tractament 1 0.55855 0.55855 1.3907 0.18817 0.2104
## Residuals 6 2.40972 0.40162 0.81183
## Total 7 2.96828 1.00000No hi ha diferències significatives entre tractaments (p-valor = 0.2104).
Espècies indicadores Abelles Franqueses
indicadores4 = multipatt(dataNMDS_abefranqueses, Metadata4$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores4)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.05
##
## Total number of species: 80
## Selected number of species: 0
## Number of species associated to 1 group: 0
##
## List of species associated to each combination:
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1No hi ha espècies indicadores d’abelles a Franqueses.
NMDS Pol·linitzadors Sabadell
dataNMDS_polsabadell <- data_shannon_MTTR %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
filter(Zona == "Sabadell") %>% dplyr::select(-Zona)
# Extreure de la taula el tractament de cada mostreig: una nova columna
Metadata5 <- dataNMDS_polsabadell %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Metadata5 <- Metadata5 %>% mutate(color=ifelse(Tractament=="NS","black","black"))# NMDS:
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)
set.seed(7)
# k=2: nombre de dimensions
NMDS_polsabadell <- metaMDS(dataNMDS_polsabadell, k=2)## Wisconsin double standardization
## Run 0 stress 0.1692459
## Run 1 stress 0.1684865
## ... New best solution
## ... Procrustes: rmse 0.1078519 max resid 0.2690531
## Run 2 stress 0.178387
## Run 3 stress 0.1720622
## Run 4 stress 0.1746865
## Run 5 stress 0.1693409
## Run 6 stress 0.1951087
## Run 7 stress 0.1783596
## Run 8 stress 0.1782982
## Run 9 stress 0.1684865
## ... Procrustes: rmse 1.220837e-05 max resid 2.176932e-05
## ... Similar to previous best
## Run 10 stress 0.1684865
## ... Procrustes: rmse 3.457522e-06 max resid 6.800015e-06
## ... Similar to previous best
## Run 11 stress 0.1684865
## ... New best solution
## ... Procrustes: rmse 2.680791e-06 max resid 5.15384e-06
## ... Similar to previous best
## Run 12 stress 0.188975
## Run 13 stress 0.1683491
## ... New best solution
## ... Procrustes: rmse 0.05383036 max resid 0.1315997
## Run 14 stress 0.1925565
## Run 15 stress 0.1693409
## Run 16 stress 0.1720621
## Run 17 stress 0.1783863
## Run 18 stress 0.1955375
## Run 19 stress 0.1900103
## Run 20 stress 0.1757075
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 17: stress ratio > sratmax
## 2: scale factor of the gradient < sfgrmincolors <- c("#74AA90", "#DB6443")
# Plot NMDS pol·linitzadors Campus
ordiplot(NMDS_polsabadell,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata5$Tractament)
for (i in seq_along(unique_groups)) {
ordihull(NMDS_polsabadell, groups = Metadata5$Tractament,
show.groups = unique_groups[i],
draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_polsabadell,display="species",col="black",air=1, cex = 1.2)
# Save the last plot
g <- recordPlot()
# Replay and save as PNG
width <- 17 # Width in inches
height <- 14 # Height in inches
res <- 500 # Resolution in dpi
# Replay and save as PNG with higher resolution and larger size
png("PolSabadell.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()## quartz_off_screen
## 2Anàlisi de similitud (ANOSIM) Pol·linitzadors Sabadell
ano = anosim(dataNMDS_polsabadell, Metadata5$Tractament, distance = "bray", permutations = 9999)
ano##
## Call:
## anosim(x = dataNMDS_polsabadell, grouping = Metadata5$Tractament, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: -0.1244
## Significance: 0.8597
##
## Permutation: free
## Number of permutations: 9999No signiticatiu (p valor = 0.8597): la composició de les comunitats no és significativament diferent entre els tractaments S i NS.
ADONIS Pol·linitzadors Sabadell
# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_polsabadell <- dataNMDS_polsabadell %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
head(Tractament_metadata_polsabadell)## Tractament
## Sabadell-Juny-2021_NS_1 NS
## Sabadell-Juny-2021_NS_4 NS
## Sabadell-Juny-2021_NS_7 NS
## Sabadell-Juny-2021_NS_8 NS
## Sabadell-Juny-2021_NS_11 NS
## Sabadell-Juny-2021_S_2 SNMDStot.dist <- vegdist(wisconsin(sqrt(dataNMDS_polsabadell)), k=2)
adonis(formula = NMDStot.dist~Tractament, data = Tractament_metadata_polsabadell, permutations = 9999)$aov.tab## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").## Permutation: free
## Number of permutations: 9999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tractament 1 0.3582 0.35823 1.026 0.09305 0.4316
## Residuals 10 3.4916 0.34916 0.90695
## Total 11 3.8498 1.00000No hi ha diferències significatives entre tractaments (p-valor = 0.4316).
Espècies indicadores Pol·linitzadors Sabadell
indicadores5 = multipatt(dataNMDS_polsabadell, Metadata5$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores5)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.05
##
## Total number of species: 137
## Selected number of species: 0
## Number of species associated to 1 group: 0
##
## List of species associated to each combination:
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1No hi ha espècies indicadores de pol·linitzadors a Sabadell
NMDS Abelles Sabadell
dataNMDS_abesabadell <- abelles_shannon_MTTR %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
filter(Zona == "Sabadell") %>% dplyr::select(-Zona)
# Extreure de la taula el tractament de cada mostreig: una nova columna
Metadata6 <- dataNMDS_abesabadell %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Metadata6 <- Metadata6 %>% mutate(color=ifelse(Tractament=="NS","black","black"))# NMDS abelles Sabadell:
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)
set.seed(7)
# k=2: nombre de dimensions
NMDS_abesabadell <- metaMDS(dataNMDS_abesabadell, k=2)## Wisconsin double standardization
## Run 0 stress 0.1316041
## Run 1 stress 0.1469069
## Run 2 stress 0.1595895
## Run 3 stress 0.1335987
## Run 4 stress 0.1657059
## Run 5 stress 0.1316041
## ... Procrustes: rmse 2.349789e-06 max resid 4.067505e-06
## ... Similar to previous best
## Run 6 stress 0.1316041
## ... Procrustes: rmse 2.165757e-05 max resid 3.932637e-05
## ... Similar to previous best
## Run 7 stress 0.1316041
## ... Procrustes: rmse 4.392255e-06 max resid 8.304797e-06
## ... Similar to previous best
## Run 8 stress 0.1336417
## Run 9 stress 0.1335987
## Run 10 stress 0.1316041
## ... Procrustes: rmse 2.786497e-05 max resid 5.137867e-05
## ... Similar to previous best
## Run 11 stress 0.1337683
## Run 12 stress 0.1316041
## ... Procrustes: rmse 9.613328e-06 max resid 1.585993e-05
## ... Similar to previous best
## Run 13 stress 0.1316041
## ... Procrustes: rmse 3.396958e-06 max resid 6.825139e-06
## ... Similar to previous best
## Run 14 stress 0.1695293
## Run 15 stress 0.1330299
## Run 16 stress 0.1335987
## Run 17 stress 0.1316041
## ... Procrustes: rmse 2.094775e-06 max resid 3.591374e-06
## ... Similar to previous best
## Run 18 stress 0.169319
## Run 19 stress 0.1673448
## Run 20 stress 0.1732531
## *** Best solution repeated 7 timescolors <- c("#74AA90", "#DB6443")
# Plot NMDS pol·linitzadors Campus
ordiplot(NMDS_abesabadell,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata6$Tractament)
for (i in seq_along(unique_groups)) {
ordihull(NMDS_abesabadell, groups = Metadata6$Tractament,
show.groups = unique_groups[i],
draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_abesabadell,display="species",col="black",air=1, cex = 1.2)
# Save the last plot
g <- recordPlot()
# Replay and save as PNG
width <- 17 # Width in inches
height <- 14 # Height in inches
res <- 500 # Resolution in dpi
# Replay and save as PNG with higher resolution and larger size
png("AbeSabadell.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()## quartz_off_screen
## 2Anàlisi de similitud (ANOSIM) Abelles Sabadell
ano = anosim(dataNMDS_abesabadell, Metadata6$Tractament, distance = "bray", permutations = 9999)
ano##
## Call:
## anosim(x = dataNMDS_abesabadell, grouping = Metadata6$Tractament, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: -0.1622
## Significance: 0.9447
##
## Permutation: free
## Number of permutations: 9999No signiticatiu (p valor = 0.9447): la composició de les comunitats no és significativament diferent entre els tractaments S i NS.
ADONIS Abelles Sabadell
# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_abesabadell <- dataNMDS_abesabadell %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
head(Tractament_metadata_abesabadell)## Tractament
## Sabadell-Juny-2021_NS_1 NS
## Sabadell-Juny-2021_NS_4 NS
## Sabadell-Juny-2021_NS_7 NS
## Sabadell-Juny-2021_NS_8 NS
## Sabadell-Juny-2021_NS_11 NS
## Sabadell-Juny-2021_S_2 SNMDStot.dist <- vegdist(wisconsin(sqrt(dataNMDS_abesabadell)), k=2)
adonis(formula = NMDStot.dist~Tractament, data = Tractament_metadata_abesabadell, permutations = 9999)$aov.tab## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").## Permutation: free
## Number of permutations: 9999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tractament 1 0.2856 0.28561 0.86626 0.07972 0.5582
## Residuals 10 3.2970 0.32970 0.92028
## Total 11 3.5826 1.00000No hi ha diferències significatives entre tractaments (p-valor = 0.5582).
Espècies indicadores Pol·linitzadors Sabadell
indicadores6 = multipatt(dataNMDS_abesabadell, Metadata6$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores6)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.05
##
## Total number of species: 80
## Selected number of species: 0
## Number of species associated to 1 group: 0
##
## List of species associated to each combination:
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1No hi ha espècies d’abelles indicadores a Sabadell.
NMDS amb presència-absència (no abundàncies) ABELLES CAMPUS
NMDS_data_abellescampus_presabs<-dataNMDS_abecampus
NMDS_data_abellescampus_presabs[dataNMDS_abecampus>0]=1# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)
set.seed(7)
# k=2: nombre de dimensions
NMDS_presabs <- metaMDS(NMDS_data_abellescampus_presabs, k=2)## Run 0 stress 0.1220077
## Run 1 stress 0.1218802
## ... New best solution
## ... Procrustes: rmse 0.05344429 max resid 0.2642867
## Run 2 stress 0.1203376
## ... New best solution
## ... Procrustes: rmse 0.07231672 max resid 0.3645329
## Run 3 stress 0.1208759
## Run 4 stress 0.119132
## ... New best solution
## ... Procrustes: rmse 0.05694998 max resid 0.3191073
## Run 5 stress 0.1191307
## ... New best solution
## ... Procrustes: rmse 0.0306343 max resid 0.1370904
## Run 6 stress 0.120841
## Run 7 stress 0.1208037
## Run 8 stress 0.1206953
## Run 9 stress 0.1207318
## Run 10 stress 0.1208075
## Run 11 stress 0.1203058
## Run 12 stress 0.120813
## Run 13 stress 0.1202204
## Run 14 stress 0.1203477
## Run 15 stress 0.1230263
## Run 16 stress 0.1232023
## Run 17 stress 0.1202727
## Run 18 stress 0.1220156
## Run 19 stress 0.1206846
## Run 20 stress 0.1212361
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 18: no. of iterations >= maxit
## 2: stress ratio > sratmaxcolors <- c("#74AA90", "#DB6443")
# Plot NMDS total
ordiplot(NMDS_presabs,type="n", xlab = "NMDS1", ylab = "NMDS2", cex.lab = 1.75, cex.axis = 1.75)
unique_groups <- unique(Metadata2$Tractament)
for (i in seq_along(unique_groups)) {
ordihull(NMDS_presabs, groups = Metadata2$Tractament,
show.groups = unique_groups[i],
draw = "polygon", col = colors[i], label = F)
}
orditorp(NMDS_presabs,display="species",col="black",air=0.01, cex = 1.2)
# Save the last plot
g <- recordPlot()
# Replay and save as PNG
width <- 17 # Width in inches
height <- 14 # Height in inches
res <- 500 # Resolution in dpi
# Replay and save as PNG with higher resolution and larger size
png("Presenciabsencia.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()Anàlisi de similitud ANOSIM presència-absència
ano = anosim(NMDS_data_abellescampus_presabs, Metadata2$Tractament, distance = "bray", permutations = 9999)
ano##
## Call:
## anosim(x = NMDS_data_abellescampus_presabs, grouping = Metadata2$Tractament, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.01252
## Significance: 0.3341
##
## Permutation: free
## Number of permutations: 9999Les composicions de les comunitats dels dos tractaments son molt similars, i el resultat no és significatiu (0.3341).
ADONIS presència-absència
# Extreure de la taula el tractament de cada mostreig: una nova columna
ADONIS_Tractament_metadata <- NMDS_data_abellescampus_presabs %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
ADONIS_Tractament_metadata## Tractament
## Campus-juliol-2020_NS_2 NS
## Campus-juliol-2020_NS_3 NS
## Campus-juliol-2020_NS_4 NS
## Campus-juliol-2020_NS_6 NS
## Campus-juliol-2020_NS_8 NS
## Campus-juliol-2020_NS_10 NS
## Campus-juliol-2020_NS_11 NS
## Campus-juliol-2020_NS_12 NS
## Campus-juliol-2020_NS_13 NS
## Campus-juliol-2020_NS_14 NS
## Campus-juliol-2020_NS_15 NS
## Campus-juliol-2020_NS_16 NS
## Campus-juliol-2020_NS_20 NS
## Campus-juliol-2020_S_1 S
## Campus-juliol-2020_S_5 S
## Campus-juliol-2020_S_7 S
## Campus-juliol-2020_S_9 S
## Campus-juliol-2020_S_17 S
## Campus-juliol-2020_S_18 S
## Campus-juliol-2020_S_19 S
## Campus-juliol-2020_S_21 S
## Campus-juliol-2020_S_22 S
## Campus-Juliol-2021_NS_2 NS
## Campus-Juliol-2021_NS_4 NS
## Campus-Juliol-2021_NS_6 NS
## Campus-Juliol-2021_NS_8 NS
## Campus-Juliol-2021_NS_10 NS
## Campus-Juliol-2021_NS_12 NS
## Campus-Juliol-2021_NS_14 NS
## Campus-Juliol-2021_NS_16 NS
## Campus-Juliol-2021_NS_18 NS
## Campus-Juliol-2021_NS_20 NS
## Campus-Juliol-2021_S_1 S
## Campus-Juliol-2021_S_3 S
## Campus-Juliol-2021_S_5 S
## Campus-Juliol-2021_S_7 S
## Campus-Juliol-2021_S_9 S
## Campus-Juliol-2021_S_11 S
## Campus-Juliol-2021_S_13 S
## Campus-Juliol-2021_S_15 S
## Campus-Juliol-2021_S_17 S
## Campus-Jun-2021_NS_2 NS
## Campus-Jun-2021_NS_4 NS
## Campus-Jun-2021_NS_6 NS
## Campus-Jun-2021_NS_8 NS
## Campus-Jun-2021_NS_10 NS
## Campus-Jun-2021_NS_12 NS
## Campus-Jun-2021_NS_14 NS
## Campus-Jun-2021_NS_16 NS
## Campus-Jun-2021_NS_18 NS
## Campus-Jun-2021_NS_20 NS
## Campus-Jun-2021_S_1 S
## Campus-Jun-2021_S_3 S
## Campus-Jun-2021_S_5 S
## Campus-Jun-2021_S_7 S
## Campus-Jun-2021_S_9 S
## Campus-Jun-2021_S_11 S
## Campus-Jun-2021_S_13 S
## Campus-Jun-2021_S_15 S
## Campus-Jun-2021_S_17 S
## Campus-Jun-2021_S_19 S
## Campus-maig-2020_NS_1 NS
## Campus-maig-2020_NS_2 NS
## Campus-maig-2020_NS_3 NS
## Campus-maig-2020_NS_4 NS
## Campus-maig-2020_NS_5 NS
## Campus-maig-2020_NS_6 NS
## Campus-maig-2020_NS_7 NS
## Campus-maig-2020_NS_8 NS
## Campus-maig-2020_NS_9 NS
## Campus-maig-2020_NS_10 NS
## Campus-maig-2020_NS_11 NS
## Campus-maig-2020_NS_12 NS
## Campus-maig-2020_NS_13 NS
## Campus-maig-2020_NS_14 NS
## Campus-maig-2020_NS_15 NS
## Campus-maig-2020_NS_16 NS
## Campus-maig-2020_NS_17 NS
## Campus-maig-2020_NS_18 NS
## Campus-maig-2020_NS_19 NSNMDSabe.dist <- vegdist(wisconsin(sqrt(NMDS_data_abellescampus_presabs)), k=2)
adonis(formula = NMDSabe.dist~Tractament, data = ADONIS_Tractament_metadata, permutations = 1000)$aov.tab## Warning: 'adonis' is deprecated.
## Use 'adonis2' instead.
## See help("Deprecated") and help("vegan-deprecated").## Permutation: free
## Number of permutations: 1000
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tractament 1 0.818 0.81797 2.1994 0.02742 0.01299 *
## Residuals 78 29.009 0.37191 0.97258
## Total 79 29.827 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El tractament explica un 2,724% de la variabilitat entre les comunitats, i el resultat és significatiu (p-valor = 0.01299).
Espècies indicadores presència-absència
indicadores = multipatt(NMDS_data_abellescampus_presabs, Metadata2$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores, alpha = 0.05)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.05
##
## Total number of species: 80
## Selected number of species: 3
## Number of species associated to 1 group: 3
##
## List of species associated to each combination:
##
## Group NS #sps. 2
## stat p.value
## Lasioglossum cf. pauxillum 0.289 0.0270 *
## Rhodanthidium septemdentatum 0.287 0.0275 *
##
## Group S #sps. 1
## stat p.value
## Lasioglossum politum 0.312 0.0109 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hi ha espècies indicadores de cada tractament!
MODELS LINEALS MIXTOS D’ANÀLISIS DESCRIPTIVES:
Abundància de Lasioglossums
# Crear matriu amb només les columnes de Lasioglossums
noms_columnes <- colnames(metadata_shannon)
genere <- "Lasioglossum"
columnes_lasioglossum <- noms_columnes[grep(paste0("^", genere), noms_columnes)]
dades_lasioglossum <- metadata_shannon[, columnes_lasioglossum]
# Agrupar totes les columnes en una sola (abundància total de Lasioglossums en cada trampa)
dades_lasioglossum$Abundancia_Total <- rowSums(dades_lasioglossum)
dades_lasioglossum_ok <- data.frame(Abundancia_Total = dades_lasioglossum$Abundancia_Total)
rownames(dades_lasioglossum_ok) <- rownames(dades_lasioglossum)
# View(dades_lasioglossum_ok)
# Extreure columnes de zona i tractament
lasioglossums_zona_tractament <- dades_lasioglossum_ok %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
mutate(Tractament = str_extract(rownames(.), "(?<=_)[A-Z]+(?=_\\d)"))
# View(lasioglossums_zona_tractament)Model lineal mixt (tractament com a variable explicativa, i zona com a variable aleatòria)
# install.packages("lme4")
# install.packages("Matrix")
model_lasioglossum <- lmer(Abundancia_Total ~ Tractament + (1|Zona), data = lasioglossums_zona_tractament)
summary(model_lasioglossum)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Abundancia_Total ~ Tractament + (1 | Zona)
## Data: lasioglossums_zona_tractament
##
## REML criterion at convergence: 642.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1440 -0.4778 -0.1408 0.1962 7.2804
##
## Random effects:
## Groups Name Variance Std.Dev.
## Zona (Intercept) 1.557 1.248
## Residual 35.226 5.935
## Number of obs: 101, groups: Zona, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.926 1.172 2.051 1.643 0.23908
## TractamentS 3.954 1.226 98.906 3.224 0.00171 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.429
Quan el tractament és NS, l’abundància de Lasioglossums és superior a 0 (Intercept: marginalment significatiu). Quan el tractament és S, s’espera que l’abundància de lasioglossums augmenti 3,954 respecte el tractament NS, i és un resultat molt significatiu (p-valor = 0.00171). Els Lasioglossums son més abundants a les zones segades (concorda amb algunes de les espècies “indicadores” de S!)
Model lineal (Tractament com a variable explicativa, sense considerar la zona)
# Convertir la variable Tractament a factor
lasioglossums_zona_tractament$Tractament <- factor(lasioglossums_zona_tractament$Tractament)
# Ajustar el model lineal
model <- lm(Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament)
summary(model)##
## Call:
## lm(formula = Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.421 -2.619 -1.421 1.381 43.579
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.619 0.754 3.473 0.000764 ***
## TractamentS 3.802 1.229 3.093 0.002576 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.985 on 99 degrees of freedom
## Multiple R-squared: 0.08811, Adjusted R-squared: 0.0789
## F-statistic: 9.566 on 1 and 99 DF, p-value: 0.002576# Crear el gráfico
plot(Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament, pch = 16)
Diversitat de Shannon d’abelles (zona com a variable aleatòria)
# Crear un dataframe de la diversitat de shannon d'abelles (calculada abans)
taula_shannon_abelles <- data.frame(Shannon = abelles_shannondiv)
rownames(taula_shannon_abelles) <- rownames(abelles_shannon_MTTR)
# Extreure dues noves columnes: tractament i zona
taula_shannon_abelles <- taula_shannon_abelles %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2])
print(taula_shannon_abelles)## Shannon Zona Tractament
## Campus-juliol-2020_NS_2 1.3862944 Campus NS
## Campus-juliol-2020_NS_3 1.6674619 Campus NS
## Campus-juliol-2020_NS_4 1.6417347 Campus NS
## Campus-juliol-2020_NS_6 1.3296613 Campus NS
## Campus-juliol-2020_NS_8 1.6674619 Campus NS
## Campus-juliol-2020_NS_10 0.7963116 Campus NS
## Campus-juliol-2020_NS_11 0.0000000 Campus NS
## Campus-juliol-2020_NS_12 1.5607104 Campus NS
## Campus-juliol-2020_NS_13 1.6434177 Campus NS
## Campus-juliol-2020_NS_14 0.7549968 Campus NS
## Campus-juliol-2020_NS_15 1.5498260 Campus NS
## Campus-juliol-2020_NS_16 0.6931472 Campus NS
## Campus-juliol-2020_NS_20 1.7478681 Campus NS
## Campus-juliol-2020_S_1 0.6931472 Campus S
## Campus-juliol-2020_S_5 1.3517840 Campus S
## Campus-juliol-2020_S_7 0.5505734 Campus S
## Campus-juliol-2020_S_9 0.0000000 Campus S
## Campus-juliol-2020_S_17 1.0986123 Campus S
## Campus-juliol-2020_S_18 1.3862944 Campus S
## Campus-juliol-2020_S_19 0.6931472 Campus S
## Campus-juliol-2020_S_21 1.1189689 Campus S
## Campus-juliol-2020_S_22 0.6108643 Campus S
## Campus-Juliol-2021_NS_2 0.0000000 Campus NS
## Campus-Juliol-2021_NS_4 0.0000000 Campus NS
## Campus-Juliol-2021_NS_6 0.5623351 Campus NS
## Campus-Juliol-2021_NS_8 1.0042425 Campus NS
## Campus-Juliol-2021_NS_10 0.0000000 Campus NS
## Campus-Juliol-2021_NS_12 1.0549202 Campus NS
## Campus-Juliol-2021_NS_14 1.0397208 Campus NS
## Campus-Juliol-2021_NS_16 0.0000000 Campus NS
## Campus-Juliol-2021_NS_18 1.3862944 Campus NS
## Campus-Juliol-2021_NS_20 0.6931472 Campus NS
## Campus-Juliol-2021_S_1 0.6931472 Campus S
## Campus-Juliol-2021_S_3 0.6365142 Campus S
## Campus-Juliol-2021_S_5 1.8662160 Campus S
## Campus-Juliol-2021_S_7 1.4369665 Campus S
## Campus-Juliol-2021_S_9 1.1973401 Campus S
## Campus-Juliol-2021_S_11 1.4388830 Campus S
## Campus-Juliol-2021_S_13 0.0000000 Campus S
## Campus-Juliol-2021_S_15 0.6931472 Campus S
## Campus-Juliol-2021_S_17 0.5091373 Campus S
## Campus-Jun-2021_NS_2 1.4735024 Campus NS
## Campus-Jun-2021_NS_4 1.0986123 Campus NS
## Campus-Jun-2021_NS_6 0.8675632 Campus NS
## Campus-Jun-2021_NS_8 2.3393717 Campus NS
## Campus-Jun-2021_NS_10 1.3862944 Campus NS
## Campus-Jun-2021_NS_12 2.0947290 Campus NS
## Campus-Jun-2021_NS_14 1.4925477 Campus NS
## Campus-Jun-2021_NS_16 1.2130076 Campus NS
## Campus-Jun-2021_NS_18 1.5607104 Campus NS
## Campus-Jun-2021_NS_20 0.0000000 Campus NS
## Campus-Jun-2021_S_1 1.6094379 Campus S
## Campus-Jun-2021_S_3 1.5498260 Campus S
## Campus-Jun-2021_S_5 1.0397208 Campus S
## Campus-Jun-2021_S_7 0.5623351 Campus S
## Campus-Jun-2021_S_9 1.5810938 Campus S
## Campus-Jun-2021_S_11 1.5595812 Campus S
## Campus-Jun-2021_S_13 1.4388830 Campus S
## Campus-Jun-2021_S_15 0.4101163 Campus S
## Campus-Jun-2021_S_17 1.3321790 Campus S
## Campus-Jun-2021_S_19 1.2424533 Campus S
## Campus-maig-2020_NS_1 0.6682485 Campus NS
## Campus-maig-2020_NS_2 0.9611277 Campus NS
## Campus-maig-2020_NS_3 0.6365142 Campus NS
## Campus-maig-2020_NS_4 1.6094379 Campus NS
## Campus-maig-2020_NS_5 0.6365142 Campus NS
## Campus-maig-2020_NS_6 1.8845210 Campus NS
## Campus-maig-2020_NS_7 0.6931472 Campus NS
## Campus-maig-2020_NS_8 1.3862944 Campus NS
## Campus-maig-2020_NS_9 1.3321790 Campus NS
## Campus-maig-2020_NS_10 0.9502705 Campus NS
## Campus-maig-2020_NS_11 1.3862944 Campus NS
## Campus-maig-2020_NS_12 0.6931472 Campus NS
## Campus-maig-2020_NS_13 1.3114313 Campus NS
## Campus-maig-2020_NS_14 1.3296613 Campus NS
## Campus-maig-2020_NS_15 1.4978661 Campus NS
## Campus-maig-2020_NS_16 0.9502705 Campus NS
## Campus-maig-2020_NS_17 0.6365142 Campus NS
## Campus-maig-2020_NS_18 1.3862944 Campus NS
## Campus-maig-2020_NS_19 0.9251291 Campus NS
## Franqueses-Maig-2021_NS_1 0.6931472 Franqueses NS
## Franqueses-Maig-2021_NS_2 0.5004024 Franqueses NS
## Franqueses-Maig-2021_NS_3 0.6931472 Franqueses NS
## Franqueses-Maig-2021_NS_4 1.6674619 Franqueses NS
## Franqueses-Maig-2021_NS_5 1.0986123 Franqueses NS
## Franqueses-Maig-2021_S_1 0.0000000 Franqueses S
## Franqueses-Maig-2021_S_2 1.0986123 Franqueses S
## Franqueses-Maig-2021_S_4 0.6931472 Franqueses S
## Sabadell-Juny-2021_NS_1 1.7478681 Sabadell NS
## Sabadell-Juny-2021_NS_4 0.0000000 Sabadell NS
## Sabadell-Juny-2021_NS_7 1.3862944 Sabadell NS
## Sabadell-Juny-2021_NS_8 1.2424533 Sabadell NS
## Sabadell-Juny-2021_NS_11 1.3862944 Sabadell NS
## Sabadell-Juny-2021_S_2 0.6365142 Sabadell S
## Sabadell-Juny-2021_S_3 0.5623351 Sabadell S
## Sabadell-Juny-2021_S_5 1.4995094 Sabadell S
## Sabadell-Juny-2021_S_6 0.6931472 Sabadell S
## Sabadell-Juny-2021_S_9 0.6931472 Sabadell S
## Sabadell-Juny-2021_S_10 0.5004024 Sabadell S
## Sabadell-Juny-2021_S_12 1.0986123 Sabadell S# Fer el model
model_shannonabelles <- lmer(Shannon ~ Tractament + (1|Zona), data = taula_shannon_abelles)## boundary (singular) fit: see help('isSingular')## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Tractament + (1 | Zona)
## Data: taula_shannon_abelles
##
## REML criterion at convergence: 162.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0268 -0.7267 0.1090 0.7752 2.3611
##
## Random effects:
## Groups Name Variance Std.Dev.
## Zona (Intercept) 6.132e-18 2.476e-09
## Residual 2.842e-01 5.331e-01
## Number of obs: 100, groups: Zona, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.08059 0.06771 98.00000 15.959 <2e-16 ***
## TractamentS -0.13912 0.10984 98.00000 -1.267 0.208
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.616
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Quan el tractament és NS, el valor de l’índex de Shannon és significativament superior a 0 (Intercept). Quan el tractament és S, s’espera que l’índex de Shannon disminueixi 0,4605 respecte el tractament NS. És un resultat significatiu (p-valor = 0,011055). Per tant, la diversitat de Shannon en abelles és significativament inferior en el tractament segat.
# Plot del model
plot(fitted(model_shannonabelles), resid(model_shannonabelles),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")
Riquesa total
# Preparar les dades
riquesa_aux <- as.data.frame(riquesa)
rownames(riquesa_aux) <- paste(riquesa_aux$Mostreig, riquesa_aux$Tractament, sep = "_")
riquesa_aux <- subset(riquesa_aux, select = -Mostreig)
# Extreure columna de la zona
taula_riquesa <- riquesa_aux %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1])
print(taula_riquesa)## Tractament riquesa Zona
## Campus-Juliol-2021_NS NS 20 Campus
## Campus-Juliol-2021_S S 21 Campus
## Campus-Jun-2021_NS NS 45 Campus
## Campus-Jun-2021_S S 34 Campus
## Campus-juliol-2020_NS NS 39 Campus
## Campus-juliol-2020_S S 24 Campus
## Campus-maig-2020_NS NS 45 Campus
## Franqueses-Maig-2021_NS NS 20 Franqueses
## Franqueses-Maig-2021_S S 11 Franqueses
## Sabadell-Juny-2021_NS NS 16 Sabadell
## Sabadell-Juny-2021_S S 12 Sabadell# Fer el model
model_riquesa <- lmer(riquesa ~ Tractament + (1|Zona), data = taula_riquesa)
summary(model_riquesa)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: riquesa ~ Tractament + (1 | Zona)
## Data: taula_riquesa
##
## REML criterion at convergence: 70.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7954 -0.4064 -0.1774 0.6667 1.1012
##
## Random effects:
## Groups Name Variance Std.Dev.
## Zona (Intercept) 86.50 9.301
## Residual 74.49 8.631
## Number of obs: 11, groups: Zona, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 26.171 6.627 2.885 3.949 0.0311 *
## TractamentS -9.501 5.237 7.267 -1.814 0.1109
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.373TractamentS: quan el tractament és S la riquesa disminueix (-9.501), però no és un resultat significatiu.
# Plot del model
plot(fitted(model_riquesa), resid(model_riquesa),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")
Riquesa PER TRAMPA d’abelles al Campus (data de mostreig com a factor aleatori)
# Preparar les dades
dades_abellescampus <- dades_abelles[dades_abelles$Zona == "Campus", ]
riquesatrampa_abellescampus <- dades_abellescampus %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Riquesa = n_distinct(ID),.groups = "drop")
riquesa_abellescampus <- as.data.frame(riquesatrampa_abellescampus)
rownames(riquesa_abellescampus) <- paste(riquesa_abellescampus$Mostreig, riquesa_abellescampus$Tractament, riquesa_abellescampus$Trampa, sep = "_")
riquesa_abellescampus <- subset(riquesa_abellescampus, select = -Mostreig)
# Extreure columna de la data de mostreig:
riquesa_abellescampus <- riquesa_abellescampus %>%
mutate(Data_Mostreig = str_extract(rownames(.), "(?<=-)[^-_]+-[0-9]+"))
head(riquesa_abellescampus)## Tractament Trampa Riquesa Data_Mostreig
## Campus-Juliol-2021_NS_2 NS 2 1 Juliol-2021
## Campus-Juliol-2021_NS_4 NS 4 1 Juliol-2021
## Campus-Juliol-2021_NS_6 NS 6 2 Juliol-2021
## Campus-Juliol-2021_NS_8 NS 8 3 Juliol-2021
## Campus-Juliol-2021_NS_10 NS 10 1 Juliol-2021
## Campus-Juliol-2021_NS_12 NS 12 3 Juliol-2021# Fer el model
model_riquesa_abellescampus <- lmer(Riquesa ~ Tractament + (1|Data_Mostreig), data = riquesa_abellescampus)
summary(model_riquesa_abellescampus)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Riquesa ~ Tractament + (1 | Data_Mostreig)
## Data: riquesa_abellescampus
##
## REML criterion at convergence: 265
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8026 -0.7445 -0.1710 0.4959 2.9434
##
## Random effects:
## Groups Name Variance Std.Dev.
## Data_Mostreig (Intercept) 0.6885 0.8297
## Residual 4.4396 2.1070
## Number of obs: 61, groups: Data_Mostreig, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.0378 0.6041 2.8500 6.684 0.00801 **
## TractamentS -0.3864 0.5427 57.1316 -0.712 0.47932
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.414En el tractament S la riquesa d’abelles per trampa disminueix, però no és significatiu.
# Plot del model
plot(fitted(model_riquesa_abellescampus), resid(model_riquesa_abellescampus),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")
Riquesa PER MOSTREIG d’abelles al Campus (data de mostreig com a factor aleatori)
# Preparar les dades
riquesa_aux <- dades_abellescampus %>%
group_by(Mostreig, Tractament) %>%
summarise(riquesa = n_distinct(ID))## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.riquesa_aux <- as.data.frame(riquesa_aux)
rownames(riquesa_aux) <- paste(riquesa_aux$Mostreig, riquesa_aux$Tractament, sep = "_")
riquesa_aux <- subset(riquesa_aux, select = -Mostreig)
# Extreure columna de la data de mostreig:
riquesa_aux <- riquesa_aux %>%
mutate(Data_Mostreig = str_extract(rownames(.), "(?<=-)[^-_]+-[0-9]+"))
riquesa_aux## Tractament riquesa Data_Mostreig
## Campus-Juliol-2021_NS NS 12 Juliol-2021
## Campus-Juliol-2021_S S 17 Juliol-2021
## Campus-Jun-2021_NS NS 25 Jun-2021
## Campus-Jun-2021_S S 19 Jun-2021
## Campus-juliol-2020_NS NS 27 juliol-2020
## Campus-juliol-2020_S S 14 juliol-2020# Fer el model
model_riquesa_abellescampus2 <- lmer(riquesa ~ Tractament + (1|Data_Mostreig), data = riquesa_aux)## boundary (singular) fit: see help('isSingular')## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: riquesa ~ Tractament + (1 | Data_Mostreig)
## Data: riquesa_aux
##
## REML criterion at convergence: 27.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5484 -0.3180 0.2212 0.5530 0.9401
##
## Random effects:
## Groups Name Variance Std.Dev.
## Data_Mostreig (Intercept) 0.00 0.000
## Residual 36.33 6.028
## Number of obs: 6, groups: Data_Mostreig, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 21.333 3.480 4.000 6.130 0.00359 **
## TractamentS -4.667 4.922 4.000 -0.948 0.39672
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.707
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')TractamentS: quan el tractament és S la riquesa disminueix (en 5.583), però NO és un resultat significatiu (p-valor = 0.24677).
# Plot del model
plot(fitted(model_riquesa_abellescampus2), resid(model_riquesa_abellescampus2),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")
Diversitat (Shannon) d’abelles al Campus (data de mostreig com a factor aleatori)
# Crear taula auxiliar: abundàncies d'abelles per ID, per cada mostreig del Campus
# Filtrar dades per Campus
dades_abellescampus <- dades_finals[dades_finals$Zona == "Campus", ]
# Crear taula de dades d'abundància només per abelles
abundanciabelles_ID2 <- dades_abellescampus %>%
mutate(Mostreig = paste(Zona, Mes, Any, sep = "-")) %>%
filter(Ordre == 'HymenopteraAbella') %>%
group_by(Mostreig, Tractament, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.# View(abundanciabelles_ID)
#Crear una taula auxiliar de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux2 <- cast(abundanciabelles_ID2, Mostreig + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux2 <- as.data.frame(taula_aux2)
taula_aux2[is.na(taula_aux2)] <- 0
#Agrupar Mostreig i tractament en una nova variable (Agrupacio)
data_shannon_abellescampus <- taula_aux2 %>%
mutate(Agrupacio = paste(Mostreig, Tractament, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament)
#Definir els rownames amb la variable Agrupacio (mostreig i tractament)
rownames(data_shannon_abellescampus) <- data_shannon_abellescampus$Agrupacio
data_shannon_abellescampus <- data_shannon_abellescampus %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# head(data_shannon_abellescampus)# Càlcul de l'índex de Shannon i preparació del dataframe
shannondiv_abellescampus <- diversity(data_shannon_abellescampus)
# Crear un dataframe de la diversitat de shannon d'abelles
taula_shannon_abellescampus <- data.frame(Shannon = shannondiv_abellescampus)
rownames(taula_shannon_abellescampus) <- rownames(data_shannon_abellescampus)
# Extreure dues noves columnes: tractament i data mostreig
taula_shannon_abellescampus <- taula_shannon_abellescampus %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>%
mutate(Data_Mostreig = str_extract(rownames(.), "(?<=-)[^-_]+-[0-9]+"))
print(taula_shannon_abellescampus)## Shannon Tractament Data_Mostreig
## Campus-juliol-2020_NS 2.611229 NS juliol-2020
## Campus-juliol-2020_S 1.905899 S juliol-2020
## Campus-Juliol-2021_NS 2.090031 NS Juliol-2021
## Campus-Juliol-2021_S 1.806531 S Juliol-2021
## Campus-Jun-2021_NS 2.701045 NS Jun-2021
## Campus-Jun-2021_S 2.080060 S Jun-2021
## Campus-maig-2020_NS 2.199104 NS maig-2020# Fer el model
model_shannonabellescampus <- lmer(Shannon ~ Tractament + (1|Data_Mostreig), data = taula_shannon_abellescampus)
summary(model_shannonabellescampus)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Tractament + (1 | Data_Mostreig)
## Data: taula_shannon_abellescampus
##
## REML criterion at convergence: 1.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.03570 -0.47827 0.03651 0.57598 0.80378
##
## Random effects:
## Groups Name Variance Std.Dev.
## Data_Mostreig (Intercept) 0.03686 0.1920
## Residual 0.02461 0.1569
## Number of obs: 7, groups: Data_Mostreig, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.4004 0.1240 3.9717 19.363 4.42e-05 ***
## TractamentS -0.5097 0.1248 2.2656 -4.083 0.0444 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.398TractamentS: quan el tractament és S, s’espera que l’índex de Shannon disminueixi 0,5097 respecte el tractament NS. És un resultat significatiu! (p-valor= 0.0444).
# Plot del model
plot(fitted(model_shannonabellescampus), resid(model_shannonabellescampus),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")
MODELS LINEALS D’ANÀLISIS FUNCIONALS
Inter-tegular span (ITS) (model lineal mixt)
Abelles totals, amb Zona com a factor aleatori
#Obrir la base de dades amb la columna de les mesures de la ITS i filtrar dades per abelles
abelles_ITS <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
abelles_ITS <- abelles_ITS %>%
mutate(data_mostreig = paste(Mes, Any, sep = "-"))
head(abelles_ITS)## # A tibble: 6 × 27
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 159 Campus Juliol 7 2021 10 NS <NA>
## # ℹ 19 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>, data_mostreig <chr># Model lineal mixt ITS (mida)
abelles_ITS_filtrades <- abelles_ITS[!is.na(abelles_ITS$Inter_tegular_span),]
model_abelles_ITS <- lmer(`Inter_tegular_span` ~ Tractament + (1|Zona) + (1 | data_mostreig), data = abelles_ITS_filtrades)
summary(model_abelles_ITS)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | Zona) + (1 | data_mostreig)
## Data: abelles_ITS_filtrades
##
## REML criterion at convergence: 1425.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3887 -0.6319 -0.1823 0.1260 5.6755
##
## Random effects:
## Groups Name Variance Std.Dev.
## data_mostreig (Intercept) 0.03817 0.1954
## Zona (Intercept) 0.01395 0.1181
## Residual 0.45922 0.6777
## Number of obs: 683, groups: data_mostreig, 6; Zona, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.69083 0.12062 1.57957 14.018 0.012 *
## TractamentS -0.40651 0.06206 592.97166 -6.551 1.24e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.230La ITS de les abelles disminueix -0.40651 quan el tractament és S (abelles més petites quan està segat). És significatiu (p-valor = 1.24e-10)
# Plot del model
plot(fitted(model_abelles_ITS), resid(model_abelles_ITS),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")
Abelles del Campus, amb data_mostreig com a factor aleatori
# Filtrar les dades només per abelles del Campus
abellescampus_ITS <- abelles_ITS[abelles_ITS$Zona == "Campus",]
#Definir nova variable: "data_mostreig", combinació de mes i any de mostreig
abellescampus_ITS <- abellescampus_ITS %>%
mutate(data_mostreig = paste(Mes, Any, sep = "-"))
head(abellescampus_ITS)## # A tibble: 6 × 27
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 159 Campus Juliol 7 2021 10 NS <NA>
## # ℹ 19 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>, data_mostreig <chr># Fer el model
model_abellescampus_ITS <- lmer(`Inter_tegular_span` ~ Tractament + (1|data_mostreig), data = abellescampus_ITS)
summary(model_abellescampus_ITS)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | data_mostreig)
## Data: abellescampus_ITS
##
## REML criterion at convergence: 1196.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3101 -0.5834 -0.2132 0.1382 5.4645
##
## Random effects:
## Groups Name Variance Std.Dev.
## data_mostreig (Intercept) 0.04629 0.2151
## Residual 0.42778 0.6540
## Number of obs: 594, groups: data_mostreig, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.54523 0.11430 3.29968 13.519 0.000529 ***
## TractamentS -0.36675 0.06534 537.59881 -5.613 3.19e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.243Molt significatiu: la ITS de les abelles (mida) disminueix -0.36675 quan el tractament és S. Les abelles son més petites quan està segat.
# Plot del model
plot(fitted(model_abellescampus_ITS), resid(model_abellescampus_ITS),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")
ITS abelles per espècies totals:
# Crear nova base de dades: només una espècie per cada tractament.
especies_ITS <- abelles_ITS %>%
distinct(ID, Tractament, .keep_all = TRUE)
especies_ITS <- especies_ITS %>%
mutate(data_mostreig = paste(Mes, Any, sep = "-"))
# Dialectus -- dialectus
head(especies_ITS)## # A tibble: 6 × 27
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 159 Campus Juliol 7 2021 10 NS <NA>
## # ℹ 19 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>, data_mostreig <chr># Fer el model (no podem considerar la zona com a factor aleatori ja que hem eliminat espècimens!)
model_especies_ITS <- lmer(`Inter_tegular_span` ~ Tractament +(1|Zona) + (1|data_mostreig), data = especies_ITS)
summary(model_especies_ITS)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | Zona) + (1 | data_mostreig)
## Data: especies_ITS
##
## REML criterion at convergence: 203.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3366 -0.6971 -0.2923 0.5604 3.1878
##
## Random effects:
## Groups Name Variance Std.Dev.
## data_mostreig (Intercept) 0.07046 0.2655
## Zona (Intercept) 0.07477 0.2734
## Residual 1.00207 1.0010
## Number of obs: 70, groups: data_mostreig, 6; Zona, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.2408 0.2593 2.1034 8.642 0.0112 *
## TractamentS -0.1192 0.2619 65.7501 -0.455 0.6505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.318No significatiu
# Plot del model
plot(fitted(model_especies_ITS), resid(model_especies_ITS),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")
ITS abelles per espècies del campus
# Crear nova base de dades: només una espècie per cada tractament.
especiescampus_ITS <- abellescampus_ITS %>%
distinct(ID, Tractament, .keep_all = TRUE)
# Dialectus -- dialectus
head(especiescampus_ITS)## # A tibble: 6 × 27
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 159 Campus Juliol 7 2021 10 NS <NA>
## # ℹ 19 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>, data_mostreig <chr># Fer el model (no podem considerar la zona com a factor aleatori ja que hem eliminat espècimens!)
model_especiescampus_ITS <- lmer(`Inter_tegular_span` ~ Tractament + (1|data_mostreig), data = especiescampus_ITS)
summary(model_especiescampus_ITS)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | data_mostreig)
## Data: especiescampus_ITS
##
## REML criterion at convergence: 168.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4350 -0.7065 -0.2774 0.5608 2.6641
##
## Random effects:
## Groups Name Variance Std.Dev.
## data_mostreig (Intercept) 0.07751 0.2784
## Residual 0.91037 0.9541
## Number of obs: 60, groups: data_mostreig, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.0006 0.2138 5.6703 9.358 0.000117 ***
## TractamentS -0.1076 0.2634 56.2365 -0.408 0.684543
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.465No significatiu.
# Plot del model
plot(fitted(model_especiescampus_ITS), resid(model_especiescampus_ITS),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")
Substrat de nidificació (GLMM binomial)
Nest_type en funció del tractament de les abelles totals (zona i data_mostreig = factor aleatori)
abelles_ITS<-abelles_ITS[!is.na(abelles_ITS$Nesting_type),]
abelles_ITS["Nius"] = 0
abelles_ITS$Nius[abelles_ITS["Nesting_type"]=="S"]="S"
abelles_ITS$Nius[abelles_ITS["Nesting_type"]!="S"]="Not_S"
# Taula de contingència
taula_contingencia1 <- table(abelles_ITS$Tractament, abelles_ITS$Nius)
print(taula_contingencia1)##
## Not_S S
## NS 71 324
## S 14 291# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia1)
ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
labs(x = "Tractament", y = "Nº", fill = "Nesting_type") +
ggtitle("Nest type totals") +
scale_fill_manual(values = c("S" = "#3D405B", "Not_S" = "#F2CC8F")) +
theme_minimal()
# Model:
abelles_ITS$Nius <- factor(abelles_ITS$Nius)
# Nomes Zona o tambe mostreig, com a factors aleatoris
model_abelles_nestt <- glmer(Nius ~ Tractament + (1|Zona) + (1|data_mostreig), data = abelles_ITS, family = binomial)## boundary (singular) fit: see help('isSingular')## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Nius ~ Tractament + (1 | Zona) + (1 | data_mostreig)
## Data: abelles_ITS
##
## AIC BIC logLik deviance df.resid
## 493.3 511.5 -242.6 485.3 696
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.7820 0.2091 0.2282 0.4715 0.5146
##
## Random effects:
## Groups Name Variance Std.Dev.
## data_mostreig (Intercept) 3.795e-02 0.1947978
## Zona (Intercept) 1.277e-10 0.0000113
## Number of obs: 700, groups: data_mostreig, 6; Zona, 3
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4843 0.1645 9.024 < 2e-16 ***
## TractamentS 1.6262 0.3450 4.714 2.43e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.442
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Significatiu! Les abelles del tractament S tenen més probabilitat (1.6262) de fer nius al sòl que de fer nius a altres llocs que no siguin el sòl. És un resultat significatiu! (p-valor = 2.43e-06).
# Plot del model
plot(fitted(model_abelles_nestt), resid(model_abelles_nestt),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")
Nest_type en funció del tractament de les abelles del Campus (data_mostreig = factor aleatori)
abellescampus_ITS<-abellescampus_ITS[!is.na(abellescampus_ITS$Nesting_type),]
abellescampus_ITS["Nius"] = 0
abellescampus_ITS$Nius[abellescampus_ITS["Nesting_type"]=="S"]="S"
abellescampus_ITS$Nius[abellescampus_ITS["Nesting_type"]!="S"]="Not_S"
# Taula de contingència
taula_contingencia2 <- table(abellescampus_ITS$Tractament, abellescampus_ITS$Nius)
print(taula_contingencia2)##
## Not_S S
## NS 63 291
## S 13 243# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia2)
ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
labs(x = "Tractament", y = "Nº", fill = "Nesting_type") +
ggtitle("Nest type Campus") +
scale_fill_manual(values = c("S" = "#3D405B", "Not_S" = "#F2CC8F")) +
theme_minimal()
# Model:
abellescampus_ITS$Nesting_type <- factor(abellescampus_ITS$Nius)
model_abellescampus_nest <- glmer(Nesting_type ~ Tractament + (1|data_mostreig), data = abellescampus_ITS, family = binomial)
summary(model_abellescampus_nest)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Nesting_type ~ Tractament + (1 | data_mostreig)
## Data: abellescampus_ITS
##
## AIC BIC logLik deviance df.resid
## 440.2 453.4 -217.1 434.2 607
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5162 0.2291 0.2360 0.4689 0.4997
##
## Random effects:
## Groups Name Variance Std.Dev.
## data_mostreig (Intercept) 0.02712 0.1647
## Number of obs: 610, groups: data_mostreig, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4968 0.1764 8.484 < 2e-16 ***
## TractamentS 1.5004 0.3860 3.887 0.000101 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.512Significatiu! Les abelles del tractament S tenen més probabilitat (1.5004) de fer nius al sòl que de fer nius a altres llocs que no siguin el sòl. És un resultat significatiu! (p-valor = 0.000101).
# Plot del model
plot(fitted(model_abellescampus_nest), resid(model_abellescampus_nest),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")
Lèctia (GLMM binomial)
Lecty en funció del tractament de les abelles totals (Zona i data_mostreig = factors aleatoris)
# Taula de contingència
taula_contingencia3 <- table(abelles_ITS$Tractament, abelles_ITS$Lecty)
print(taula_contingencia3)##
## Oligolectic Polylectic
## NS 142 249
## S 38 263# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia3)
ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
labs(x = "Tractament", y = "Nº", fill = "Lecty") +
ggtitle("Lecty totals") +
scale_fill_manual(values = c("Polylectic" = "#3D405B", "Oligolectic" = "#F2CC8F")) +
theme_minimal()
# Model:
abelles_ITS$Lecty <- factor(abelles_ITS$Lecty)
model_abelles_lecty <- glmer(Lecty ~ Tractament + (1|Zona) + (1|data_mostreig), data = abelles_ITS, family = binomial)## boundary (singular) fit: see help('isSingular')## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Lecty ~ Tractament + (1 | Zona) + (1 | data_mostreig)
## Data: abelles_ITS
##
## AIC BIC logLik deviance df.resid
## 589.1 607.3 -290.6 581.1 688
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1853 -0.8175 0.1735 0.5793 1.2233
##
## Random effects:
## Groups Name Variance Std.Dev.
## data_mostreig (Intercept) 4.057 2.014
## Zona (Intercept) 0.000 0.000
## Number of obs: 692, groups: data_mostreig, 6; Zona, 3
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.5679 0.8603 1.823 0.0684 .
## TractamentS 0.6398 0.2749 2.328 0.0199 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.119
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Les abelles tenen una major probabilitat de ser polilèctiques (0,6398) quan estan sotmeses al tractament S. És un resultat significatiu (p-valor= 0,0199), la qual cosa suggereix un efecte del tractament en la lecticitat de les abelles.
# Plot del model
plot(fitted(model_abelles_lecty), resid(model_abelles_lecty),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")
Lecty en funció del tractament de les abelles del campus (data_mostreig = factor aleatori)
# Taula de contingència
taula_contingencia4 <- table(abellescampus_ITS$Tractament, abellescampus_ITS$Lecty)
print(taula_contingencia4)##
## Oligolectic Polylectic
## NS 127 223
## S 20 232# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia4)
ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
labs(x = "Tractament", y = "Nº", fill = "Lecty") +
ggtitle("Lecty Campus") +
scale_fill_manual(values = c("Polylectic" = "#3D405B", "Oligolectic" = "#F2CC8F")) +
theme_minimal()
# Model:
abellescampus_ITS$Lecty <- factor(abellescampus_ITS$Lecty)
model_abellescampus_lecty <- glmer(Lecty ~ Tractament + (1|data_mostreig), data = abellescampus_ITS, family = binomial)
summary(model_abellescampus_lecty)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Lecty ~ Tractament + (1 | data_mostreig)
## Data: abellescampus_ITS
##
## AIC BIC logLik deviance df.resid
## 468.8 482.0 -231.4 462.8 599
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2618 0.0518 0.1706 0.5796 1.2245
##
## Random effects:
## Groups Name Variance Std.Dev.
## data_mostreig (Intercept) 6.291 2.508
## Number of obs: 602, groups: data_mostreig, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.2170 1.3182 1.682 0.0926 .
## TractamentS 0.6375 0.3271 1.949 0.0513 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.091Les abelles tenen una major probabilitat de ser polilèctiques quan estan sotmeses al tractament S. És un resultat marginalment significatiu (p-valor = 0.0513), la qual cosa suggereix un cert efecte del tractament en la lecticitat de les abelles.
# Plot del model
plot(fitted(model_abellescampus_lecty), resid(model_abellescampus_lecty),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")